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Chapter 1 
Introduction 


1.1 Background 

The transonic flight regime has proven to be very difficult to analyze because of the 
inherent nonlinearities associated with the equations governing the flow at transonic 
speeds, even in their most simplified form. It is well known that, due to these 
nonlineanties, small changes in geometry or freestream conditions do not necessarily 
produce small changes in the flow solution around the vehicle. Particularly for 
unsteady problems, which is the case in any flutter analysis, the shock motions 
and variations in the shock strength can be significant even for small oscillations of 
the body. In some cases shock-boundary layer interaction phenomenon can strongly 
influence the point of flow separation. Besides being another nonlinear effect, it also 

means that the position of separation is moving as the shock location and strength 
vaxy with the body motion. 

These combined problems have, over the years, meant that no methods were 
developed which can treat general transonic flows. On the other hand, despite all 
its difficulties, the transonic regime is very important in aeronautical applications. 
Methods that would be able to perform routine aeroelaatic analyses are necessary, 
smce current transport aircraft usually cruise in this regime and military aircraft 
maneuver in it. Even missiles, or launch vehicles in general, that are designed to 
acquire a very high speed in a reasonably short time, have to pass safely through 
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transonic speed,. Such launch vehicles are the focal subject of the p^ent invest!- 

gation. . . 

From the point of view of the aeroelasticUn the n«d to take mto cons.derat.on 

the inherent nonlinearities of transonic flow is specify annoying because it prevent, 
hi. accustomed u* of superpositioa Of courae, this is no, to say that one can never 
find useful linearised solutions in the transonic caw. For instance, the work o an- 
dahl HI show, that if, the reduced frequencim are Ugh enough, linefeed equarions 
do a good job of predicting the unsteady aerodynamic force, on streamlin wmg, 
and bodies. Another example of linearized tmnwnie solutions be found map- 
plication, of the iniiM method (e. g. , Ballhau. and Goorjian ® and Nixon ) . 
This method is a more limited form of superposition in which unsteady «rodynmmc 
solution, are given - Unear perturbation, about nonlimmr steady state solut.ona 
For smaU ampUtude. of body motion and small shock morions, tins method 
proved to be an acceptable approach to aeroelastic analysis. 

For the problem, treated here, however, no linearisation is po».bl. smee the non- 
iitmar nUure of the aerodynmnic forces is e-entml to capture the physical phenom- 
ena involved. This work undertakes to analyze the aerod-tic st.bd.ty of bdhsric 
veUdes during their tnu-onic ph»* of flight. Although the methodology devel- 
oped can be used for any flight vehicle, the reader will recognize that the p^tculm 
problems being addressed are most likdy to occur on bdlistic vehicle, (boosters), 
especiaUy those carrying so-caUed hammer heed payloads. In such cases, the pro 
lem can be compounded by the existence of relatively Urge region, of separated 


Since the meaning of a hammerhead payhmd may not be familUr to dl realms, 
it is important to mention that throughout this work the term will denote all thorn 
Uunch vehicle, where the payload has a larger diameter than the adjacent booster 
stage. It is clear that such a configuration will always be aviated wtth some form 
of boatt.il right at the front of the booster, which create, the pomibUity that at 
leut some of the vehicle could be immersed in a repo- of separated flow. 

The interest in such veUcle configurations, as weU as in blunt geometries in 
general, is reUted to the fact that they are sometimes very attractive options for 
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the designer. Aeroelastic problems have, however, been observed in the past on 
certain vehicles with these configurations t 4 ’ 5 ' 6 ) when passing through the transonic 
regime. A study of nonlinear hammerhead effects is presented in considerable detail 
by Woods ud Ericsson ^ for an Atlas-Able IV launch vehicle, which was one of 
the configurations that actually experienced flutter problems during transonic flight 
in the early 1960’s. Another example where the hammerhead could have been a 
possible source of aeroelastic instability was the original Seasat-A launch vehicle, 
as described by Ericsson and Reding W. In this case, the solution adopted com- 
pletely eliminated the hammerhead by enclosing both the payload and the Agena 
upper stage in a fairing with a diameter equal to the first-stage Atlas booster di- 
ameter. Besides true flutter cases, some failures in early launch vehicles may have 

been associated with buffeting during the transonic phase of flight, as docribed bv 
Rainey W . 

The classical approach to the analysis of the aeroelastic phenomena present in 
the transonic phase of flight of ballistic launch vehicles has been to use experi- 
mental data for the unsteady aerodynamic pressures together with some simplified 

structural-dynamic reputation of the vehicle. It is typified by the pfa 

preyed in Reference, [4], [5] and [6] . Actually, not only in the analysis of 
launch vehicle problem, but in any situation where tramonic and separated flow, 
are present, the usual practice has been to rely on empirical, or semi-empirkaL 
methods for the treatment of the aerodymunic terms, as can be seen from Refer- 
ence, [7] through [101 • The main difficulty with these analy** is exactly that they 
require the use of experimental data for the unsteady aerodynamic pressures in or- 
er to formulate the aeroelastic problem. In some cases, since unsteady pressure 

not usudly available, sternly data me adapted in a quasi-steady fimhion. 

Whole of rely ’ ng antirely on experimental aerodynamic data may be 
highly undesirable in many situations. For instance, in the emly stags, of design of 
a new vehicle, or when some modification, to an masting vehicle have to be made 
to accommodate a payload with a larger diameter than the booster’s uppermost 
stage, such experimental data are not available. The problem is usually 
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building wind tunnel models, testing them, and modifying the design if aeroelas- 
tic problems appear. This can be very expensive during the design cycle, not to 
mention that it is perhaps another reason the aeroelastician is usually seen as a 
policeman as mentioned by Ashley and co-authors l 11 *. Thus, the motivation for 
the present work includes the development of the capability for calculating, as part 
of the solution, the aerodynamic environment the vehicle is subjected to. Given the 
necessary aerodynamic tools, the aeroelastic stability analysis can be performed. 
Eventually, this computational procedure will become sufficiently cost-efficient to 
be incorporated in the design cycle, at least in the transonic regime where the usual 

linearized methods break down. 


1.2 Computational Approach 

The recent progress in the field of Computational Fluid Dynamics (CFD) has al- 
lowed the simulation of transonic flowfields through the use of finite difference, or 
finite volume, techniques. However, much more current research is focussed on 
steady transonic calculations than on unsteady ones. The reason for this is that 
unsteady calculations, and their aeroelastic applications, require that the equations 
be solved in a time- accurate mann er. As a result, the time-step size? that can be 
taken in the time integration of the equations are severely restricted and increase 
considerably the computer cycles required for solution. 

There is also the question of at which level of approximation to perform those 
computations. For some problems, it is sufficient to consider the transonic small 
disturbance equation, whereas in others a Navier-Stokes formulation would be ap- 
propriate. Interesting surveys in this regard, indicating also the stage of develop- 
ment at the time of publication, are presented by Peterson ^ and by McCroekey, 
Kutler and Bridgeman ^^1. It is no surprise that most of the work that can be found 
in the literature is restricted to airfoil or wing flows. The aerodynamic formulation 
is restricted to the transonic small disturbance or to the full potential equations, 
since these formulations are less computationally d e mand i ng than more complex 
ones based on the Euler or the Navier-Stokes equations. 
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Son* examples of applications of CFD methods to unsteady transonic aerody- 
namic calculations and to aeroelastic analyses, where the aerodynamic formulation 
is based on the small disturbance or the full potential equations, can be found in 
References [14]-[30]. The computational requirements, even in those cases, are by 
no means trivial. Certainly they are already within acceptable limits of present com- 
puters, as the work in these references shows. The study of transonic aileron buzz 
by Steger and Bailey ^ can be considered one of the “damicaT example, of the 
use of CFD for aeroelastic analysis: even before 1980 it employed a Navier-Stokes 
formulation for the Sow solver. A few other cares, where similar complex formula- 
irons are used, can also be found in the literature. For example, Reference, [32] and 
[33] include applications of Euler equations, but still to airfoil or wing problems. 

The main difficulty when one starts to consider transonic aeroelastic problems of 
bod.es is that usually the disturbance, are large enough that a potential formulation 
" ”° k>n *' r fWtkennore, the description of the aerodynamic phenomena 

rmportant to the aeroelastic analysis, such as the topology of flow separation or 
s oc - undary layer interactions, may be beyond the scope of potential methods 
In other words, one ha. to nesort to the Euler or Navier-Stokes equation, in order to 
appropriately simulate the physical flow feature, involved W. For these cue, very 
trie has actually been published, although the idea of how one should proceed to 
perform aeroeUstic analysis is somewhat well established The computational 
requrrements, both in terms of time and storage, are the major challenge, since 
®trd systems that would support a Navier-Stokes solution for typical vehicle 

shape, are bound to be hrrge, and CPU time, associated with unsteady tnumonio 
Navier-Stokes solutions are also substantial. 


1.3 The Present Method 

The approach followed in the present work uses CFD techniques to perform aeroe- 
lastrc analyses of launch vehicle configuration, by coupling the structural-dynamic 

equations representing the vehicle with an unsteady flow solver appropriate for the 
physical situation being treated. 
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Th« mechanism underlying the aeroeUstic instabilities observed in flight involves 
■ • of phase lags between the aerodynamic forces and the motion. In sue 

circumstances, these tad. can do positive work on the emulating vehicle. For Uunch 
vehicle configurations, the accurate calculation of the* lag. involve, rmolving the 

shock-boundary layer motion., mid shock-induced flow *pmabon. The 
. ,. ■ particularly important for hammerhead configurations, where 

Of flow separation » P“*‘ ciD have a definite influence on 

z 2= tSL = i « — - — • ~ ; 

vehicles being emphasized here, tPe correct u -m 

‘ h % N "C o7“«'Z sets of equation, is formed by integrating both of 

themsimultaneously in time mid ensuring that the data gmmrated by ooeset - «d 

in the next time step of the other. The aerodynmnic equation. *££££> 
a t > n< 4 their solution provides the forcing terms for the aeroelas 

’ZZZ Which typically me weighted hkph of •>- pr«*uredistn^ution alongthe 

body surface. Solution of the structural-dynamic equations gives ne 

shape of the body, and so the boundary condition, for the eerodynamic solution 

th ' Tte XHl «I of the aerodynamic equations de*rv«s a momde- 
uuled discussion. A Motion would be considered rime-eccuie m a numerical sens 
ifaUUhe scale, supported by the computation* m»h are being accuracy 
which in practical term, implies tha, the CFL numb* is at meet of the mder ne 
everywhme in the computation* domain. Since the CFL, or Courant number can 
- the ratio of the Imigth sc*« obt*n«i by the product of *me 
characteristic velocity time, the time step over the length scale determi ned by the 
SI it is clear that this is a very restrictive condition in term, of the mmomurn 
step, sizes that could be u*d- f nis is not th. meaning we intend to assign to 
expremion in the present work. Here, the Motion is being called t.m«-.ccur. 
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in an aerodynamic sense, meaning that all the important aerodynamic phenom- 
ena of interest are being resolved accurately. Since the time scales associated with 
the structural-dynamic phenomena we will be studying here are much larger than 
those associated with the pure aerodynamic phenomena, as we will demonstrate 
later, by enforcing time-accuracy in an aerodynamic sense we automatically have 
time-accuracy in an aeroelastic sense. 

Since full elution of the Navier-Stolrm equations around complex geometri „ 

“ Stm beyond our '°®putation»I capabilities, the method employed here solve, 
the thin layer approximation to the Reynold,. Avenged N.vier-Stolce. equations, 
where turbulence closure is obtained by using an algebraic eddy viscosity model 
Before the aeroelastic analysis could be started, the aerodynmnic equation, must 
be marched m time in order to obtain an initud stead, state solution for the flow 
around the body, assuming that it doe. not deform under the lemds. Once this 
mrtml solution* obtained, on. can start to oteallate the body and, from then on, 
t e aerodynamic solution must be time-accurate since the flow solver will be coupled 
o the structural-dynamic equation, ss previously dmcribed. By tracing the growth 

or decay of. perturbed oscillation, the aeroelmtic stability of a given configuration 
can be ascertained. 


In the following chapters, the formulation of the problem and the detailed de- 
scnption of the method will be presented, together with the result, obtained in 
e vanous c^ analysed. The theoretical formulation that underlie, the present 
approach wdl be presented in Chapters 2 and 3 , Chapters 4 through 6 will describe 
several applications of the method, both for unsteady aerodynamic calculation, ss 

Ch * Pter 7 possible contri- 

of the present research, present its conclusions, and digues some mcom- 
mendations for future work. 
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Chapter 2 


Formulation of the Aerodynamic 
Problem 


2.1 The Navier-Stokes Equations 

In order to correctly represent the physical phenomena present in this problem, the 
appropriate set of flow equations .to be used are the Navier-Stokes equations. These 
equations can be written in differential form as : 

continuity equation: 




conservation of momentum equations: 


energy equation: 


pEE.-vp+v.n-^ 

= % + V. [f.J-?l + p<5.u + Q 

Dt ot 


( 2 - 1 ) 


( 2 . 2 ) 


(2.3) 


where the symbol D/Dt indicates the substantial or material derivative, the vector 

Q represents the body forces, and Q is the heat addition. 

The set of equations above requires some constitutive relations in order to form 
a closed system of equations. Starting with the equation of state, a perfect gas 
is assumed such that one can write 


p = pRT — (7 - 1) p«. 


(2.4) 
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where the specific internal energy of the fluid, e„ is given by 

'■ = C - T (2.5) 

The heat flux vector is given by the Fourier law of host conduction 

?--/c(W) (2.6) 

The components of the viscous stress tensor, ? , can be obtained from 

_ (du{ duj\ 2 

T - 1=l ‘{dT, + ai:)-3^ V - ir > S -‘ (2.7) 

where we are considering an isotropic fluid mid also that the coefficient of bulk 

™coe.ty A « simply given by: A = -ip . Finally, the total enthalpy of the fluid 
i 2 , is denned as 

Bmh + j|u| J we, + £+I|s|» (2.8) 

where for a perfect gas one can also write the enthalpy h~CT 

Equation, 2.i, 2.2 and 2.3 are written in the non<ons«rvntion form This may 
create numerical problems when computing flow quantities acme, a shock wav, 
smce the us. of the non-conservation form can cause low, or creation, of mass and 

"trr" ^ ^ fr0 “ 1 ■" PO ”' ^ very important 

dZTw h ^ form or divergent form Before 

however, we will introduce the assumptions that in the problems being 

^ “ n ° forces (<3 = o) mid no heat addition (<J = 0) It 

‘‘I" 0 " 4 lh “ ^ we me refeling to her, 

field effect, such a. grav.ty, f or example. The introduction of the* assumptions 

I ,ttP “ th “ P ° ml “ * h ' However, since they would 

Tf- ^ <hAt “* "“° n f “ —• *^v, assumptions^ physic^' , 

“ iodldt» b * m * lre4led ^ C ° ntribU,i0n ° f ^ *“ 

J“* iD ““ n ° toti “’ ““ repeated-index implying 

ion, and remembering that the substantial derivative is given by 


n* 37" + | 


(2.9) 
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one obtains the following equations : 



dp dp dxij _ 

a + “'a ^ +, ’^' 0 

(2.10) 


[a* ami __ap . a** 
p [ dt u, a*>] a*,- 

(2.11) 


-It* -'3 

(2.12) 

Defining the total energy per unit of volume, e , as 



t s P [e. + Jtiju;] 

(2.13) 

the equations 

can be written in the conservation-law form in 

the following way, 

continuity: 

Qt dxj Ky *' 

(2.14) 

momentum: 

d(pui) d(pUiUj) dP__?ZiL = o 

dt dxj dxi dij 

(2.15) 


en * rSr gj+gj-[(e + ri»y- T «“‘ + «l ” 0 (216) 

Full solution of Equations 2d4-2.16 Mound complex gsometrie. is still beyond 
our computational capabilities, since the number of grid points required to capture 
d the scales of a turbulent flow at flight Reynold, numbm would be Pr»b.b.tnr« 
The alternative approach consist, then in averaging the flow variable, such that 
they can be expremed M a mean flow quantity plus a sen^memi perturbatron. In 
particular, if a m. weighted averaging is u«d, it is prwible to rewrite the ><«•*» 
in term, of the averaged quantity Ernest in the same form - Equatron. 2.14-2. . 

A generic quantity x is averaged such that 


where: 7 = Jr f£+ Tp z dt 

Here the period of integration T, should be small compared to the rime 


scale for 
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variations in the mean flow quantities but should be large enough to provide some 
averaging over the high frequency turbulent scales. 

The flow variables can be written as: 

P = P + / 

pV% = f!>U t + (pUif 

e = 7+e' 

p = (2.18) 

r 0 = ^ + r'- 

p# = 

ph = pK+( P hy 

However, if averaged velocities ere nectary, for example, eome other form of av- 
«agmg le required, since the quantity that was averaged above is the momentum 
So, for quantities like velocity and enthalpy, Mother average is M 

u< s Will (2.19) 

h = JKj-p 

and with these new definitions one could write 


= U. + U? 

h = h + h" 


(2.20) 


It IS important to note, however, that perturbation quantities such as u" and h" do 
not necessarily have zero mean value. 

If we perform the averaging process on Equations 2.14-2.16 , the Mowing gov- 
erning equations are obtained, 
continuity: 


d , 

dt + fc j(P5>) = ° (2.21) 


ll 
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momentum: 


energy: 


+£<>«*>+ Sr 


( 2 . 22 ) 


3 


d 

dxj 


5i (r y -^) + <^-¥ j] 


3xj 


= 0 


(2.23) 


One can observe that the new terms are: 
pu'-u'j called Reynolds stress terms; 

-^K called Reynolds heat flux terms; 

,, ( r _ called Reynolds dissipation terms. 

V? 3 J- that these new terms would require some additional 

Tt U evident at this point that tnese new «»u« 

T . , however, that as one tries to derive a new equation 

closure equations. It so happens, however, , , 

Closure en ir if another new term is introduced, 

¥ calculate say, the second order tensor pv t u } , 

to calculat , y, constitutes the so called hoWence 

r ly - — “ •«— * - ° f u ~ 

Z“ZL, S, leas, another new tensor one order Usher than the previous one is 

"utrnate approach, which we shall follow here, is the Boumines, concept of 
. [38, 391 h by the turbulent mixing is modelled by upgrading 

t:ZZ:lr vi^colien, hy son, entity unmlly cUled 
eiscosity coefficient. We will postpone the diwossion of the P«t.eular turbulence 
model being used her. until a later action. By now, it is sufficient to say ha th 
turbulent mixing will be modelled by an eddy visccmty 

viscoeity coefficient that appear, on the definition of the v«ous strms term, wrll 
be formed as « 4 \ 

H„e p, is th. molecule viscosity coefficient and p, is the eddy visco-ity coUBcient. 

Similarly, the coefficient of thermal conductivity « can be obtamed as 

_ C,p, CjP, (2.25) 

k <- Kt + «« - -pT* + p rt 
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where Pr is the Prandtl number and Pr, is the turbulent Prandtl number 

With the introduction of the model shove mentioned into Equation. 2 21-2 23 
they can now be expressed only in terms of the averaged qualities, i. , all thJ 
terms that involved Perturbation, (primed and double-primed terms) arc in some 
sens, replaced by the model. From now on, w, shall drop all the bar, and titda, in 
he vanaUes, jus, to simplify the notation, but one should understand that w. are 

1^1 >V ' ra *' d ,U “ titi “' T1 “ ' , '“ ,iOT5 «“ th “ b « «•»*•« the 

averaged variables as 

continuity: 



ar + a^(p“i) = ° 

(2.26) 

momentum: 

d d 

dt + dx ~ ( ' pUiUj + P 6 ** “ r *>) = 0 

(2.27) 

energy: 

de d 



dt + di “ l(c + u i ~ r *3 u . + qj] = 0 

(2.28) 

where 




(2.29) 



(2.30) 


e * = ~ " \ u i u j (2.31) 

The above set of equation, constitute, the socallcd Reynolds-Averaged Navier- 

“ the - - -u z 

nares As wTJr ^ 7“ " " ti “' n lboV ' for 1 «“•—» set d coordi- 

attent' " y ' “ n ° l “ ,ir * 15 ' d " irabl '. and so we will tum our 

attentron now to the problem of coordinate transformation. 
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2.2 Coordinate Transformation 

Before proceeding with our discussion, it is useful to rewrite the flow governing 
equations in the matrix form they me usually seeninCFD applications. So, Equa- 
tion, 2.26-2.28 can be rewritten in strong conservation-law form, stiU in cartesian 

coordinates, as: 0Q dE 9G _ Q (2.32) 

where the vector of conserved quantities Q is 

P 

pu 
pv 

pw 
e 

and the flux vectors E , F and G are 

pu 

pU 2 + P “ Trr 
pUV — T*, 
pUW — T gt 

{e + p-T„)u - TgyV - T t ,w + qr 

pv 

pUV — Tgg 

PV 2 + P~T n 

pv w — r v , 

(e + p - r„) v - r„u - r yt w + q„ 
pw 

pUW — Tgg 
pVW — T,, 
pW 2 + p— Ttt 

(e + p - r M ) w - r XM u - r yM v + q, 







14 



CHAPTER 2. FORMULATION OF THE AERODYNAMIC PROBLEM 


It is convenient, for the purpose of actually implementing the flow solver code, to 
have the governing equations transformed to a general body-conforming curvilinear 
coordinate system. This will make the formulation in the code independent of 
the details of the actual topology being solved, besides making it very convenient 
for applying the boundary conditions and implementing turbulence models. The 
transformation is usually known only numerically, which does not pose any difficulty 
since what one really wants to know are the metrics of the transformation and its 
Jacobian. Following the usual procedure in the CFD literature (see, for instance, 
References[41H45]), one can convert the Navier-Stokes equations from cartesian 
coordinates to general curvilinear coordinates by means of the transformation : 


T = t 

t = {(x,y t z,t) 

T} = »7(*,y,r,<) (2.37) 

C = C(x,y,z,t) 


Using chain rule expansions, the derivatives in terms of the cartesian variables 
can be expressed in terms of the curvilinear derivatives, in matrix form, as 


a 

Tt 

9 

& 


1 *« 7* <* 

0 e. 7. Cr 

0 7, C, 

0 7. Cm 


(2.38) 


or conversely, 


* 1 •» ► * 1 1 b I 
J I _ 0 *<!/<*< I £ I 

* (2 - 39) 

£ ' 0 x < y< *< J l h J 

Due to the special form of the transformation matrices above, the Jacobian of the 
transformation, J , can be simply written 


J- I* (*,7.0 /*(*,*,*)! 
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or in the form it is actually used in computation 

J~ l = |d(x,y, *)/£({,»?, C)l 


(2.40) 


Substituting the transformation expressions given by Eq. 2.38 mto Eq. 2.32 , one 
can obtain: 


dQ dQ dQ^ r dQ dE dE_ dE_ 

dF dF .dF .dG dG_ . dG_ _ Q 

+ c, ac +£ 'af c * 


(2.41) 


Multiplying the above equation by }-' , and completing the « .n> .uch that at 

least some of them can be written in divergence form, we have 

L ( r'Q ) + ^ [J- 1 {(,Q + i.E + US + <* G >1 

+. L [r* M + n.E + n,F + -i.G)] + ^ i-r' «■« + <- £ + C,f + C.G)] 

(242) 

+ E + ac ' c *l 

+° [If + ^ + ac } = 0 

If we follow the work on Ref. [ 41 J , and use Eqs. 2.3S and 2.39 , an expression 
for the Jacobian of the transformation can be found as 

J = (x < y,z c + x„y<z< + x<y<*„ - *<y<*n “ *nVt*< ~ x cV»i z <) 

Similarly, the foUowing metric relations can be found 

z, = J (y,z c - y c*„) e, * J ( x < x n - x n x <) 

e, - J (*„y< - *<Vn) J (**< - v<*<) 

rjy = J (x<z< - x<z<) rj, = J (x ( y< - *<y<) 

C* - J (y<*, “ Vf> z i) (v = J ~ X ^ 

Cm-J (*<yn - x »>y«) 


(2.43) 
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~ Xr t* ~ y^y “ z rii (2.44) 

* 7 * “ VrVy ““ z rT)t 

Ct = ” X rCx ~ J/rC* "• ^rCi 

Working with the above metric relations and the Jacobian expression, one can prove 

that all the terms inside the curly brackets in Eq. 2.42 are identically zero. If we 
define, then, 


$ = J~'Q = J~ x 


P 

pu 

pv 

pw 

e 


* = J~ l (ttQ + {*E + t v F + S M G) 
T 5 j~ l (mQ + rj,E + rj y F + T!,G) 
& - J ~'(<ttQ+<;sE+<;yF+( t G) 


(2.45) 


the governing equations can be rewritten, still in strong conservation-law form, but 
^ ® general set of curvilinear coordinates, as: 


5F dT dU 

dr + aT + ^ + ac =0 (2.46) 

This matrix equation, or its equivalent set of scalar equations, is essentially what 
is implemented for the flow governing equations on the code developed. The above 
still does not really explain the details of the numerical implementation of these 
equations, and a later section will be dedicated to address this issue. 


2.3 Turbulence Model 

As previously mentioned, the concept of an effective viscosity coefficient is used 
here in order to model the turbulent mixing by upgrading the molecular viscosity 
coefficient by the so called eddy viscosity coefficient. In the present work, the eddy 
viscosity coefficient, p t , is obtained from the two-layer Baldwin and Lomax t 46 l 
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algebraic model, which is implemented here in the usiml -ay for waU-bounded 
shear layers !47,48 'as 

= ( (fit) inner • ^ '/crossover (2.47) 

/i ‘ \ outer . »/ > '/crossover 

Here n is the rareilmeer distance normal from the wall and 1 cro«over is the snmllrat 
value of , at which the inner and outer formulations give equal value.. 

The inner region uses the Prandtl-Van Driest formulation, such that 

Winner = P 1 ' M 

„here | u | is the magnitude of the local vorticitjr vector, and the length scale 1 is 
obtained as [ / » ? + M 

I = **/ [l _cxp \ A+ j] 

and i 7 + = (rfy/p^z) /fi*,k = 04 ’ A ” 26 . clauser formulation and is 

The formulation for the outer region is similar to a Clauser 

PV “ ^ Wouter = KC^hFkuM (2 '“ 9) 

where / r. 

F^...h. smaller of 

The values of and are obtained from the function 

p(i) = im[ 1 -«*p(- a*)] 

where F™ is the maximum value of F(q) in the profile, and w is the value of 
q at which it occurs. The function F K ,M . <*Ued the Klebanoff intemnttency 

factor, is given by ,-i 

«•*“&)] 

Finally, C*/ is the difference between maximum and minimum total velocity mag- 
nitude. in the profile (note that the minimum total velocity b sere . for boundary 
layers), and the constants used have the following values : K - 0.0168 , * ’ 

C^k = 0.25 , and Ckm — 0.3 . 


) 


18 


CHAPTER 2. FORMULATION OF THE AERODYNAMIC PROBLEM 

2.4 Numerical Implementation 

It is always good practice to perform some form of nondimensionalizatioa on the 
equations being used, even if it merely is to ensure more generality in the formu- 
lation. The choice of dimensionless parameters is somewhat arbitrary, so long as 
one does things consistently. In the present work, the density p is scaled by the 
freestream density , the cartesian velocity components ti, v and w are nondimen- 
sionalized with respect to the freestream speed of sound o w , and the total energy 
per unit of volume e is referenced to Poo<*L • 

For simplicity of notation we will not use different symbols to denote the nondi- 
mensional variables, but we will continue to use the same nomenclature as before. 
The reader should keep in mind, however, that we will be refering to the nondimen- 
sionalized variables for the remainder of this development. Equation 2.46 can still 
be written as 

#7 3T dT BU 

dr + d( + dr) + ft ;** 0 

where the vector of conserved quantities <7 is still written as previously defined. 
The reader should remember, though, that we are now referring to the nondimen- 
sional quantities, despite keeping the same nomenclature. It is instructive at this 
point, however, to rewrite the flux vectors with all their components in the gen- 
eral curvilinear coordinate system and to include the modifications caused by the 
nondimensionalization process. Thus, they can be written as 

pU 

puU + pt:-*fc. (r lx 4 x + + r„0 

pvU + rf, - + r y ,i M ) > (2.50) 

P*>U + pt 2 -*fc (r„e, + + r„is) 

(c+p)U - p£ t - fa (0'Z, + + /3,£ m ) i 
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pV 

puV + pq f — (TxsI* + TgyTJy + T xzHm) 

T = J~ l ^ pvV + prjy — (t**1* + T n1* T *»1*) 

pwV + PH j — ( r **1* + T **1* + r **1*) 

(e+p)V -pqt-*jg- (0*1* + 0*1* + 0**7*) J 

pW 

puw + pc* - ( r *»c* + r *»c» + r **c*) 

v = j * 1 { <*>w + pc* - ( r *»c* + r »»^ + r »*&) 

pwW + pQt — (TxmC* + r »*Cir + r *»C*) 

. (e+p)W-p(,-‘fc(l>JU + 0,C, + M ) 

In the above equation., i» the freeetream Mach number defined a. 

M„=- 

a™ 


(2.51) 


(2.52) 


(2.53) 


where M» - ia the magnitude of the freeetream velocity vector. He 

is the Reynold, ntimber, given in the usual way by 

^ Poof^oo^O 

Rt 


(2.54) 


Moo 


where p* is the freestream (laminar) viscosity coefficient and i» the reference 
length. It should also be mentioned that any viscosity coefficient that appears 
on the formulation is nondimensionalized with respect to Poo- The contravanant 
velocity components, U, V and W, are defined as 


U - & + + Z* w 

V - It + 1* u + 1* v + 1* w 
W = C* + C«tt + C»» + C*w 

The p x , P y and 0, terms are given by 

Q x = T„U + TzyV + r xs w - q* 

0y = r^u + + r yJ u? - q, 

0, = r x ,u + r yt v + r„w - q t 


(2.55) 


(2.56) 
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All the other variables that appear in the above equations have been previously 
defined, and again we call the reader’s attention for the fact the dimensional vari- 
ables should be just replaced by their nondimensional counterparts in those def- 
initions. It is interesting to note that the three dimensionless numbers usually 
encountered when nondimensionalizing the Navier-Stokes equations, namely Mach 
number, Reynolds number and Prandtl number also appear. Although not explic- 
itly shown in the equations immediately above, the Prandtl number arises in our 

case through the way in which we have defined the components of the heat flux 
vector, Qy &nd . 

The Beam and Warming implicit approximate factorization scheme^ 49, 5 °] is 
used for the solution of the finite difference equations. The spatial derivatives in 
Equation 2.46 are approximated using three-point, second-order central differencing, 

and the implicit Euler method is used for the time march. With this method, we 
can write 

(2 . 57) 

where by O (At 3 ) we mean that this is a first order method in time, and the su- 
perscript n indicates at which instant of time the quantity should be evaluated. If 
Equation 2.46 is substituted in the above, we obtain 

7 VH =3'-Af(H+^ + ||) + O (Al 1 ) (2.58) 

Smce the flux vectors are nonlinear functions of the vector of conserved quantities, 
m order to solve for ^ n+1 , while maintaining the order of accuracy of the method 

the nonlinearity is removed by a local Taylor series expansion about . This’ 
process yields 

- E" + ^(5" , -5’) + o( a* 1 ) 

r+1 = r*+B”(<r +, -<r)+o(A*’) (2.59) 

zr* 1 = zr + c- - ij") + o (ai“) 

where A, B and C are the Jacobian matrices given by 

. du 
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B = 


C = 


dT 

dG 

a5 


Substituting the linearized expressions into Eq. 2.58 one obtains 

d 




(2.60) 


where I is the identity matrix. 

In order to mike the solution algorithm more coat effective, an approximate 
factorization of the three^imeoaional operator into three ^dimensional operators 
U introduced. Following the work of Beam and Warming , and Pulliam , a 
term of the form 


M 


dB n dA n . dC n dA* , dB n dC n 

+ 9 - <s> T 


drj d( + ac di drj d( 


\ . A dB" 9C n dA n ] ( 

ac * J V 


At 


can be added to the left-hand side of unfactored equation. Furthermore, one should 
note that the added term is of order At 3 and therefore does not alter the order of 
accuracy of the method. With the above term, the equation can be factored such 
that the algorithm can be written in the so-called delta form as 


(2.61) 


( J * A % B ') (' + u k c ') (' * ai & A ’) ^ " 

where A, is a forward difference operator in time, such that 

A, O' = TT*' - 3" 

With the introduction the three-point second order central differencing to ap- 
proximate the space derivatives, the final form of the algorithm can be written 
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a a 

L*L(L( = R{ + R, + R( (2.62) 

where the above operators are defined as 

L ( - (j+Att'A'-erAtJ-'VtAfJ-AtMooRe-'Zj-'Mfj) 

L„ = (I + At6 v B* - - AtM'.Re-'^J-'Mfj) 

L < = {l + A<$<C n - e/At.7 -1 V<A< J - AtA/ 00 Re' 1 ^‘j~ 1 M < n j) (2.63) 

= -AtStlT - eeAtJ- 1 ( V ( A<) 3 
= -Atf n F*-« £ A<7- l (V,A,) 3 J$" 

= -A^ZT-^AtJ-^VcAJ 3 

Here ^ , 6„ and are central difference operators; V ? , V, and V { are backward 
difference operators; and A< , A, and A ( are forward difference operators in the 
£" > T 7" i Mid C-directions, respectively. For example, 

S &*> ~ 2 [^* +1 

The 5; and are midpoint operators used to maintain a compact three point 
second order accurate central difference scheme when differencing the viscous terms 

m the left-hand side. As previously described, the A, is a forward difference operator 
in time. 

The Jacobian matrices were split such that A, B and £ contain, respectively, 
the inviscid terms of A, B and C, whereas M (t M, and M c contain the viscous 
terms. Expressions for these matrices are given in Appendix A. A few extra words 
may be important to clarify the need for midpoint operators when working on 
the viscous terms. The left-hand side matrices, as defined in Equations 2.63, are 
block tridiagonal matrices where each block is 5 x 5. However, the viscous terms 
themselves, i. e. , the components of the M matrices, already involve derivatives of 
the velocity components in the case of the viscous stress terms and of the internal 
energy in the case of the heat flux terms. Therefore, if we insist in using central 
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differences, we need to use midpoint operators in order to avoid a five point stencil, 
which would cause the operators to become full-block pentadiagonal. With the 
midpoint operators we can keep the matrices block tridiagonal while retaining the 
order of accuracy of the method. 

Although linear stability analysis shows that the fully implicit algorithm is un- 
conditionally stable, stability bounds are encountered in practice. The problem 
comes about because of the nonlinear interactions in the convection terms of the 
momentum equations. The problem can be understood easily if one thinks in terms 
of waves interacting. When two waves interact, one of the products is a new wave 
of higher frequency, which frequency is the sum of the original ones. This frequency 
cascading will, at some point, exceed the resolution capacity of a finite mesh, with 
the result that these frequencies either alias back into the lower frequency range or 
pile up at the high frequency aide. If not controlled, this process may cause serious 
inaccuracies and quite possibly numerical instability. The usual way the problem 
is hwnd W is to introduce some form of numerical dissipation into the algorithm, 
with an error level that should not interfere with the accuracy of any viscous effects 
being captured in the solution. Some numerical schemes, genetically known as up- 
wind schemes, intrinsically have this numerical dissipation built into the scheme by 
the way they use one sided differences. Central difference schemes, however, do not 
have this n ume rical dissipation built into the scheme, and so it must be explicitly 
added in order to control the nonlinear instabilities above described. 

The procedure adopted here consists of introducing a constant -coefficient, fourth- 
order artificial dissipation in the right-hand side operators, and constant-coefficient, 
second-order artificial dissipation in the left-hand side operators. In the right-hand 
side, or explicit side, the amount of artificial disspation introduced is controlled by 
the coefficient , and on the left-hand side, or implicit side, it is controlled by 
the coefficient «/ (see Equations 2.63). Ideally, one would like to use fourth-order 
artificial dissipation on both implicit and explicit sides, because the use of different 
orders of numerical dissipation schemes does introduce a small numerical error on 
the solution. However, once again the issue of solution efficiency dictates that we 
settle for the second-order numerical dissipation on the left-hand side, because a 
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fourth-order scheme would cause a five point stendl and so block pentadiagonal 
matrices. 

On. possible way to dal with the problem of having different order, of artificial 
dissrpation schema on the implicit and explicit side, is to use the soiled i.ago- 
oaf algorithm, described by Pulliam and Chatxwee I 43 ! . The idea of the diagonal 
algorithm is based on the diagonaliration of the inviscid left-hand side matrices as 
described by Warming, Beam and Hyett [51) . The*- authors show that the Jaco- 
btan matrices A, B and C have real eigenvalues and a complete set of eigenvectors, 
and so can be diagonalized. It is clear that for viscous calculation, the algorithm 
cannot be rigorously applied, since in order to diagonalize the left-hand side oper- 
ator, we need to neglect the implicit viscous terms. However, result, by P ulliam 
and Sieger show that for steady viscous flows mid in convection domi-.-l un- 
steady flow, the diagonal algorithm, a. described above, produce, very good results 
Furthermore, it .flow, the us. of the fourth-order artificfal dfa.ip.tion scheme in 
the left-hand srde operator, because, despite the fact that this makes the matrices 
block pentadiagonal, now each block is a five by five diagonal matrix, making the 
inversion process fairly inexpensive even for a block pentadiagonal system. In the 
p»«nt work some attempt Was mad. to program the diagonal algorithm, however 
there were question. concerning the accuracy of the algorithm in the cam of self. 

exertrf unsteady calculation. * , which is exactly what our aerodastic analyses 
&rc. Hence the idea wu &b&odoned. 

, ^'’rMl 11 n0t ' XpKciU) ' mrat ‘ on ' d “ the Previous equations, frees tream sub- 

freefon -s perfonned in the flux vectors when computing the right-hand side 

enn. m Equation 2.42 . The reason for this comes is that, when arbitrary curvi- 

near coordinate, and general finite difference, are rued, there are small 

errore mtrodu«d in the calculation of the metric, of the transformation which may 
cause the code to be unable to reproduce the freestream (or an uniform flow). By 
perfomung the freestream subtmetion in the flux vectors, we ensure the capability 
of recovering freestream and reduce the overall emx of the method 

Since a vfacou. formulation is used here, at body wall. w. have noslip bound- 
ary conditions, ■■ '•,« = « = u, - 0 for steady problems, or « - *„ , „ » y , 

Pulliam, T.H., personal communication, May 1986. 
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&n dw = z r for unsteady cases. Here *„ Vr and z T are the cartesian components 
of the grid velocity due to body motion. We assume also adiabatic walla mid zero 
normal pressure gradient at body walls. Uniform freestream is enforced at the up- 
stream and fir-field lateral boundaries. The treatment given to the downstream 
boundary depends on the value of the freestream Mach number. For a supersonic 
freestream, the downstream boundary is extrapolated from interior values. For a 
subsonic freestream, according to characteristic relations, there are four character*- 
tic propagating downstremn mid one propagating upstream. To be consistent with 
that pressure is then fixed at the downstream boundary (at ite freestream «lue) 
and the other qualities are mctr.pol.ted from interior vjues. From a stnctly nu- 
merical point of view, it should be mentioned that all the boundary condition, are 
treated explicitly in the premrnt work. This mem* that they me applied using the 

information available at the present time step. 

There are some imports point, concerning the boundary conditions that should 
be further discua«d. In an attempt to save computational points in the gnd, the 
approach primarily used in the present work door not consider a body base, but 
rather stop, the computation at some point Jong the cylindricJ afterbody section. 
This may cause some concern when computing subsosiic Bows, since we will be as- 
suming freestream pressure at this computational exit plane whereas the flow will 
probably not be completely back to a freestream condition due to the presence of 
the body. Essentially, there is an error being introduced, but care is being exercised 
in order to make the afterbody cylindrical section long enough that whatever errors 
are introduce at the downstream boundary will not propagate upstream to the 
point of influencing the region of internet in the solution. Later, some rasults will 
be presented that illustrate this point and confirm our statement that the regions 
of interest are not being contaminated by whatever errors are introduced at the 

downstream boundary. 

Another point of concern is the treatment of the far-field lateral boundaries. 
Since we have a body at angle of attack, lift will be general* and so there is a ques- 
tion about the conservation of circulation in the far-field. To be precise, far-field 
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lateral boundary conditions should have some way of enforcing this conservation- 
of-circuiation boundary condition. However, since we are dealing with three dimen- 
sional bodies of revolution in the present work, we do not feel that the added ac- 
curacy is worth the complication of the boundary conditions. Moreover, the result, 
obtained seem to corroborate our sumptions, indicating that assuming freestream 

conditions in the upstream and far-field lateral boundaries b good enough for the 
present case. 

Fmally, the reasoning behind the zero normal pressure gradient condition at the 
wall should be addressed. The cod. b implemented with the option that n b 
the normal direction, and so this boundary condition b enforced by setting dp/d, to 
aero at the wall. Strictly speaking, for a vbcou, formulation there b no theoretical 
retiring that would lead to the zero normal pmsure gradient condition without 
the introduction of any simplifying assumption. In particular, for viscous flow, over 
curved walls tins condition is not completely comet. However, if we study mult, 
from boundary layer theory l*J mid evidence from experimental measurnments on 

tlHhT "a ° b ”' rTe llUU ““ PrOBUre " *PP rCTdm »b-ly constant throughout 
he thickness of the bounds byer. In our case, „e can abo argue thaf thb 

wJ r r *T Pti ° n U ^ *1“ Point, that are nearest to the 

w^foo.h« words, although there b a physical boundary byer, the zrn normal 
p sure condition is only being used up to the first grid surfoce off the body mid not 
broughou, the .bole byer. I, should be pointed out that Am are cases where thb 
condrtum can be ngorously demonstrated. For ,h. general cm in a N.vier-Stokcs 
formulation, hoawver, the justification for it. me b bmd <a experimental as well 
as numerical, evidence m conjunction with the insight provided by boundary byer 

For emrh time step Equation 2.62 b solved by forming the right-hand side mid 

T?" 17 ' ““ " ,h ' M,U “ Ce «f the side 

^atorn through the use of an t-U decompoeition algorithm for block tridiagonri 

27' th T Wbm “ “°” *"“** — *•'- circumferential d.rec- 

* nd “ C0 “Pl«Wy -rapped around the body for the full 360*, the 
gn is sai to be period, c in the circumferential direction. For periodic grids, it b 
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necessary to program periodic tridiagond solvere, because in this ante the lrit-hand 
side matrix is not strictly tridiagonal but has extra blocks m the upper ngh an 
the lower left comer.. Appendix B ha. more details about the handbag of penod.c 

tridiagonal matrices. .ntaimise. the amount of 

The database is structured in a pencil format wmcn mimnu 

data that has to be stored in main memory at a given time. penci 

structure concept, as well as the procedure ft* solving Equation 2.62 whem the 

ZZZ ha. sueha format, am dmcribed in detdl by Deiwer, and Rot Wd 

A lorial idea of a general body configuration in phyricsl space, together wrth ho. 

this is mapped into computation* space and showing the various block, botmdanm, 

is ormentedb. Figure 2.1 . Pencil, of data are formed by stonng s^uentmlly the 

blocks in the given sweep direction while keeping the pencil t«e restncted to one 

block dimension in the other coordinates. Incidentally, tins* a g pom 

mention that the code primarily used (a. the flow solver) * tlu. work 

a steady state version of the ARC3D code that the author, of Reference had 

,o run on a CDC Cyber 205 computer. The cod. is highly vector, 

1 makes extensive urn of random, or ssynchmnous, I/O * order to unprov, t«s 

efficiency for a problem that is definitely not core contauwd. 

The computational mmh for all the cases considered in tin. work was genera 
using algebraic methods, since for the regular geometne. <yp.calof l»unc ve c 
configurations them method, are able to provide sufficiently good gnds andthey 
are very simple to urn. It is important to realise that, since the goal ts toperf rm 
aeroelaatic analyms, the body will be deforming » the solutton proweds. 
means that some form of grid reshaping is necew-ry in order to ««-»** 
deformation. In the pr-ent approach, the complete grid - regenerated every 
to, step (when performing aeroebmtic calculation.), winch » another 
it i. important to urn Jgebmic grid generation method, winch are very f»tjf on. 
decide, to arudyxe more complex g~metri«s, »y a vehicle b °“ ’ 

probably a more elaborated grid generation scheme should be adopted. 

Both exponential and hyperbolic-tangent' 56 ' grid stretching techm^s am v* 
where necemary to clust M grid point, in the mgions of higher flow gmdrents. One 
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(a) Physical space. 



(b) Computational space. 


b<>dy CODfi * uratl °" “‘u.tntfm* the pencil (end block) dot. 
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obvious instance where such clustering must occur is dose to the body in the normal 
direction, since we need . very refinri grid in order to capture the viscous effect 

Other situation, where clustering is necemary may include Wh ^ 

m exD4c ted to occur, or region, of rapid expansion, for instance, where the body 

rcgird grid clustering is that, despite the foci the vi*ous term. 

»U three coordinate direction, in the pr-ent formulation, the appro«h . w* 
should still be comidemd within the scope of the thin layer approxunario^ t fo the 
Reynolds-Averaged Navier-Stokes equations^ . This is becai-e the gnds used me 
too coarse in the longitudind mid circumfemntid direction, to capture comp « 
the viscous effects in these directions. 
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Aeroelastic Formulation 


3.1 Equations of Motion 


The structural-dynamic equations for the vehicle are develop* by considering free- 
e. flexural vrbratron of an elongated beam with variable properties. They are cast 
m motW form. « is burned that the angle of attach renmins few enough that no 
^eral forces wrU appear, and so bending hr only one plane ha to be coLderT 
h^formulatton developed is quite general, even allowing for rigid body degree, 

f freedom, although in the application, presented here the rigid body degre7 of 
freedom were assumed to be constrained. X f 

th Tb , e “‘“‘T* ° f m ° ti0n f ° r * ® mOTl b “" “ vibration can be found in 

extensive literature (see, for instance. Reference, 1 5!) H 62] I n t—,, , 

ncenn, beam theory «veral degree, of approximation exist from whi^tl 

:rr: ti t 0 ^ ^ *■*■ # u o °* , 0 

fV ! n i ° f ““ “° n ° f th ~ — - win sim- 

to m ™ ' * h * P '* “ d “ lural vibration are available 

or from 7 *“ ° btaiMd dther structure-dynamic tat, 

we will "T a “ t * * l * meDt struct “* 1 of the vehicle. Moreover, 

mulat t ° f b ° W S ‘ ruc ‘ ural d< “P in « “ introduc* into the for- 

mulation but srmp.y take i, a a given , entity, obtain* either from expedient 

or some theoretical model. Finally, it should be clear that, although 
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aeroelutic formulation is nonlinear (as we hope to demonstrate throughout this de- 
velopment) , the strictly structural-dynamic formulation is linear and the principle 
of supcrpoeition can be applied. 

For the sake o £ completeness, end to clarify the dimension of the qualities 
involved here, the equation of motion for * beem in flexural vibration ie pwentrf. 
When the effects of shear and rotary inertia are neglected. 


&6 

m ( x ) + 



F,(x,t) 


where 4(l , t ) is the total lateral deflection of the centerline of the beam, m(x) is the 
mass per unit of length, EI(x) is the bending stiffoera, and F,(x,t) is the lateral 
external load applied to the beam (units of force per unit of length). Again, , it is 
very important to stress that the equation above is merely an example of wha a 
very simplified equation of motion for lateral vibrationof an elongated beam might 
look like. A. mentioned by Bisplinghoff and Ashley' 60 !, -.hear, rotmy inertia, and 
several other effects are automatically accounted for” when a modal superposition 

solution technique is used. , [59] 

The total deflection 6(x,t) at any station along the body can be expressed , 

uwing a modal approach, as 


$(*,<) = £«<(*) *<(*) 

iml 


(3.1) 


where * (t) are the generalised normal coordinates and * (s) are the normal modes. 
To be consistent with the aerodynamic formulation, we will assume that «(*, t) and 
,.(,) have been nondimensionalixed with respect to the same reference length f„ 
^d for the aerodynamic equations. It should be mentioned that the u« of a 
modal approach assumes that the structure and the frequencies of interest ere such 
that disturbances are felt almost instantaneously throughout the structure. In other 
words, problems with travelling structural waves, which are common in very large 
sprme structures, me of no concern for the vehicle, conmdered here. It is also clem 
that in actual application only a finite number of modes is considered. In other 
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words Eq. 3.1 should be more properly written as 

im 1 

where N represents the number of modes employed in the analysis. 

Since the normal coordinates are not coupled either elastically or inertially, the 
equation of motion for the i-th cooordinate can be written as 


m ' [*»'(*) + 2 CA?.(<) + = P-(f) 


(3.2) 

Her. the quantity have already been nondimensionalired, and the dot, indicate 

envattvee with respect to (nondimensional) time. The .mndimensional natural 
frequencies, uTj, are defined by 

(3.3) 




ence length lit ^ “ the Mh mod '> “ the refer- 

nee length used for nondnnenaionaliration of the aerodynamic equations, and «. is 

' 7 ***** 0f 3 ° und - K » “*«■«"« *o point out that the nondimensional 

equency above is not the same as the socalled reduced JrejneW 59 ' 40 -* 3 ) which 

or unsteady aerodynamic applnmtious. 
ver, since that reduced frequency, k, is usuaUy defined as 

. _ vt 0 

they are related by 

O* = M m ki 

The nondimensional gener^sed masse,, I*. can be obtain* from the actual vehicle 
mass distribution, m(x), from 


__ t m(x) 


(3.4) 

wh«e p. » the freestream density and f, is the (nondimowonal) body length. If 
«( mensional) generalized masses rfn are already known, say, from experimental 
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results, it is dear that the integration described above is not oeamqr, and the 
nondimenaionatlization is simply obtained by 

(3.5) 


m = 


TTli 


P oo^O 


It is important to emphasize that in the pr««n>t appro** the natural n , 

the structural damping coefficient. «.), the norm* mode, md the « dm nbu 
« considered input data They are turned to be known ether from tests 

other theoretical analysis. ..men. 

The real difference in the present approach, when compared to cUsstctd aeroe- 

lastic analysis, is associated with the generalised aerodynamic forces, W). " '<= 
in this case are calculated from 

P,(t) = J*' e(x,t)<t>i{x)dx ( 3 - 6 ^ 

s, ,) i, the nondimensional running norm* force acting on the vehicle, ob- 
frJtm suitable c.rcumferent.al integration, of the body pressure dts.nbuuom 

Note that the dm.nsi.nsl counterpart to «) would have units of force per urn 
ol length, and also that since w. are dealing with small angle of attack case, <(r , ) 
is numerically very close to the running lift on the body. The important pom 
that the true nonlinem character of the pre«nt analyst is somewhat MJcn . - 
genemlized aerodynamic forces, became the pre-ur. dmtnbu.ton mound^he body 
u obtained from the numeric* solution of the Navier-Stokes equattons. Not. that 
P M is identified with mode i, but h- contribution, of all nrodes through • 

In this fashion *1 trarmoruc aerodynamic nonlimmritie. me captumd by the method, 
whereas the structural-dynamic formulation is still kept very simple. 

Some point, regarding the -sumption that the angle of attack 
should be further cUrifiri. First, and probably the most important mmon why 
assumption is introduced, is the fact that the kindofbooe.cn .. are * 

in her. are prevented from attaining large angle, of a.t«k because the hmd. «lm 
structure would be tremendous. Most such vehicle, have -me mmomum angle 
of attack that, if exceeded, cause the launch to be aborted by the de.truct.on 
of the vehicle. For the Saturn V , this maximum tolerable angle of attack was 
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approximately 4 degrees* . Second, there is nothing to indicate that the l«nd of 
instability we are concerned with is a high-angle-of-attack phenomenon. Quite to 
the contrary, since the major mechanism that drives the instability is associated 
with lags in the aerodynamic forces due to the shock motion ^ 64 ^, its effect would 
be more pronounced at very s ma l l angles of attack where the stream wise location 
of the shock can actually be forward in the leeward side (when compared to the 
windward side) for a portion of the motion and backward for the rest of it. Another 
important observation is that there is no conceptual difficulty in extending the 
present method to treat cases where asymmetric separation, or yaw *ng1» exists, 
such that bending in both planes and twisting of the vehicle are to be considered. 
It is just a matter of coding the necessary geometric considerations and introducing 
the structural modes associated with the added degrees of freedom. As far as the 
aerodynamic formulation is concerned, it is valid for any angle of attack; its more 
stringent limitation is associated with the turbulence model which originally was 
derived for attached or mildly separated flows. 


3.2 Solution of Aeroelastic Equations 

Since the calculation of the flow solution at each time step, and so the evaluation 
of P^t) , is much more time consuming then the solution of the structural-dynamic 
equations, some constraints are imposed on the numerical method that can be used 
for the time integration of Equation 3.2. The idea of transforming the second 
order equation into a first order system as usually done in control theory may «**m 
attractive, but it is not very practical in this case. For instance, any implicit method 
or any predictor-corrector sequence, which would typically involve the evaluation 
of the farcing term at an instant of time ahead of the current time, would become 
prohibitive as far as computational time is concerned with present computers. On 
the other hand, the explicit Euler method can be shown to be always unstable in 
such situations by a simple linear stability analysis. 

The procedure selected to advance Equation 3.2 in time consists of a straight 
finite differencing of both the first and second time derivatives, as also done by 

Ashley, H., personal communication, July 1987. 
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Stegcr and Bailey t 31 ', Gurus wamy and Yangt 65 !, and Borland and Rizzetta f66] • 
Using second order accurate formulas, we can approximate those derivatives as 




3 


+ ?t,n— 1 

(AO 5 


9i 


2A t 


Here the subscript i refers to the mode considered and the subscript n refers to 
the time level at which the quantity is evaluated. With the above formulas, the 
expression for ?., n+ i , i. e. , the generalized deflection of mode * at tune n + 1 , is 


given by 

[2 - (At) 3 o?] ft, n - (1 - AtCA) fr-" 
tt.n+1 = (1 + At(jZ7j) — — 


(3.7) 


where, for simplicity of notation, we have denoted P,(t) =» P,(t)/W,- 

Linear stability analysis of the above scheme shows that it is conditionally stable. 
However, this really poses no constraint in the time step size for aeroelastic anal- 
ysis, because the values of At required for the stability, and more importantly the 
time-accuracy, of the aerodynamic equations alone is much smaller than whatever 
restrictions could be imposed by the conditional stability of the above scheme (for 
values of natural frequency that would be of any concern for aeroelastic stability). 
More on the stability properties of this Equation 3.7 scheme will be discussed in 
the next section. The analysis also shows that for C. - 0 the scheme is numerically 
non-dissipative, which is of course the expected result since central differences are 
being used. This is important to avoid hiding physically unstable solutions because 

of numerical dissipation introduced by the method. 

The coupling of the two sets of equations, the Navier-Stokes equations governing 
the flow behavior and the above described structural-dynamic equations governing 
the Vehicle behavior, is performed in the following way. At each time 

step, the aerodynamic equations are solved in a time-accurate fashion, and their 
solution provides the forcing terms for the aeroelastic analysis. These are weighted 
integrals of the pressure distribution along the body surface, as one can see tr-m 
Equation 3.6 . Solution of the structural-dynamic equations gives the new deformed 
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shape of the body, as well as the body oscillatm* velocities, and therefore the 
boundary conditions for the aerodynamic solution at the next time step. By tracing 

the growth or decay of a perturbed oscillation, the aeroelastic stability of a given 
configuration can be ascertained. 

Note that once the response in terms of the generalized coordinates, 9 ,’s, is 
calculated, the total deflection of the body centerline, *(*,<), can be determined. 
Since we are assuming that body cross-sections do not deform, this means that the 
total deflection of the body surface is known, and the grid can be regenerated. The 
solution process, as implemented here, recalculates the whole computational mesh 
at every time step to account for the deformation and motion of the body. This 

was done to provide greater generality to the method and enable it to treat not so 
small deflections. 

Finally, it should be stressed that, in the present implementation of the method, 
the rigid body degrees of fleedom were assumed to be somehow constrained, and 
only the elastic ones were considered. There is no additional conceptual difficulty 
m including the ngid body modes in the current code. However, results by Woods 
and Encsmn seem to indicate that the inclusion of them mode, is not critical for 
determining the aeroelastic stability of a given vehicle, for the kinds of configurations 
and phenomena dealt with in the present work. 


3.3 Stability Considerations 


la this section we intend to discuss further the numerical stability of the scheme 
implemented for the solution of the aeroelastic equations. The main objective here 
is to substantiate the claims made in the previous section regarding the stability 
•eeumey of the algorithm. The techniques used in this analysis are the same 
ones usually employed in the analysis of computational fluid dynamics schemes and 
are dmeribed in detail by Lomax*. A linear stability analysis of the scheme 
implemented for the solution of the aerodynamic equations, i. e. , the compressible 
Navmr-Stokes equations, can be found in the literature (see for instance Refer- 
ence ) and will not be discussed here. 


- ^ Umarica ^ MflthPtia in Fluid Mechanics Notes for course AA214A, 
eronautics and Astronautics, Stanford University, Autumn Quarter, 1983. 
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The numerical algorithm bring implemented here ia given by Equation 3 . 7 . Smce 
the generalized aerodynamic force are a forcing term in this <mse, for stability anal- 
yai. we should take ft,. = 0 . Conridering where e .. the sometunes 

crilrf amplification factor, and replacing this into Equation 3 . 7 , one obtarna 

(2 - At 3 o?) qi-m - (l - AtC.c r 0 <y ~ 1 g«> (3.8) 

"(1 + A*CA) 

Cancelling out the common term and solving for <r we obtam 

(2 - At 3 Q?) ± ~ (1 ~ CO 

2(1 + AlCiSTg) 


a — 


( 3 . 9 ) 


The stability requirement is given by 


W<1 


( 310 ) 


or, in words, the stability region to described by the circle of unit radius in the 

compiex^e p^te q{ ^ lbove expression for d is Probably beyond the potnt 

w. want to make here. It to interring to not. that there am 
parameters that have to be considered, namely ATO and (,. Smce for most pr«t. 
application, the structuml damping coefficient to a snmll qumrtity, rt .. mstruettve 
,o study the limiting c«e when <, to zero. In this c«« the express, on for the 

amplification factor becomes 

( 3 . 11 ) 


„ „ I f (2 - 4Po?) ± - <] 


which show, that the scheme to numriically stable foe < 2 

zero structuml dmnping. Since the structuml frequencies that would *«*»««* 
ore usuaily below 50 to 80 Hr, and also due to the fact that the 
frequencies are bring refried to the speed of sound, it » sab to state tha . vnU, 

at it. largest, be a number around unity. However, this dm- that the c<mdit,ond 
stability of the scheme petoe. no additionri constraint on the tune step, that co 
possibly be taken, bemuse the mmdmum time step required to ensure tmie-accumcy 
of the aerodynamic equation, don. is at least one order of magnitude, mid quite 
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possibly two orders of magnitude, smaller than the one determined by the above 
condition. Moreover, the expression above shows that the modulus of a is exactly 
equal to 1 for all values of < 2, which means that the numerical scheme is 
also non-dissipative v . 

In support of previous statements, we will consider now the other possibility 
that has been suggested before, i. e. , to transform the second order equation into 
a system of first order ones and then use the explicit Euler method for the time 
march. We will show that this produces an unstable numerical algorithm. Rewriting 
Equation 3.2 as a first order system produces 


f qi(t) = 

\ i'i(t) = P,(t) - U*q,(t) - 20ZJ.I /,( 


(3.12) 


where we are essentially defining a new variable s g t , an d the dots indicate 
derivatives with respect to nondimensional time, as before. Using the explicit E tiler 
method, we can write 


I j, a tori'. m „ 

1 « * Suajp* . />(,„ - - 2(.a in , 


which then produces the following set of finite difference equations 


i 9«\n+l — ?i,n + 

*Wl * (1 - 2A*CA) + 




For stability analysis, the forcing term is set to zero and we assume that qi >n +i = &qi,* 
and i/, tM>1 = <rui tn . Writing the resulting equation in matrix form, we obtain 


(» - 1) l r i _ 

-a®? <» - 1 + say*,) J \ j 


The characteristic equation is obtained by setting the determinant of the coefficient 
matrix to zero, which gives 

(<r - 1) (<r — 1 + 2A*Ci3.) - = 0 


* Lomax., H., Numerical Methods in Fluid Mechanics. Notes for course AA214A, 
Dept, of Aeronautics and Astronautics, Stanford University, Autumn Quarter, 1983. 
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If we consider again the limiting case of zero structural damping, the roots of the 
above equation can be easily obtained as 

<r u * 1 ± £*Ui 

which shows the method is always unstable for any value of parameter . 

Numerical experiments support this theoretical result, indicating that in this case 
the instability of the explicit Euler method is a very practical problem. 
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Study of Hemisphere-Cylinder 
Cases 


4.1 Introduction 

The aerodynamic code being used in the present investigation evolved from a version 
of the ARC3D code^^ that Deiwert and Rothmund had running on a CDC 
Cyber 205 computer. The code had been optimized for the Cyber 205 architecture, 
but the formulation therein did not allow for general unsteady problems since it 
was missing the metric terms associated with time derivatives. The first major 
programming task undertaken here consisted of the introduction of the unsteady 
terms into the aerodynamic equations. Before attempting any aeroelastic analysis or 
even such complex configurations as the hammerhead, it was important to check the 
new code on some problem with a simpler geometry that would allow for validation 
of the modifications introduced. 

Due to the simplicity of the geometry, while still keeping the general launch 
vehicle shape, a hemisphere-cylinder configuration was chosen to test the code. 
Such a configuration is also attractive from the standpoint that experimental or 
other computational results are available from the literature. The following sections 
will describe the cases studied with such configuration, which include steady state 
results for freest ream Mach numbers of 0.5, 0.6, and 1.5 (all cases at zero angle of 
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attack), as well as unsteady results for a rigid body pitching oscillation. It should 
be noted that for these calculations the turbulence model was turned off, so that 
they are laminar computational results. Turbulent flow cases were run only for 
the hammerhead payload configurations, which will be discussed in the subsequent 
chapters. 

4.2 Grid Generation 

The body conforming computational mesh was generated using algebraic methods, 
and the same mesh was used for all hemisphere-cylinder cases analyzed. Grid lines 
run in the longitudinal, normal and circumferential directions, and 50, 40 and 20 
grid points were used, respectively, in these directions. It should be pointed out 
that this is a fair ly coarse grid system, since again the point of these computations 
was mainly to pinpoint any possible problems with the code before attempting any 
aeroelastic solutions on more complex geometries. 

A general three dimensional view of body and grid can be seen in Figure 4.1, and 
details of typical grid planes, i. e. , longitudinal and crossflow planes, can be seen 
in Figure 4.2 . Mesh points in the normal direction are clustered near the body in 
order to capture viscous effects, and a 25% exponential grid stretching is used in this 
direction. Over the hemispherical part of the body, grid lines in the longitudinal 
direction are placed at equal angular increments, and over the cylindrical part of it, 
these lines are equally spaced. The 50 points used in the longitud in al direction are 
distributed such that 15 of them are over the hemispherical part of the body, and the 
other 35 are located along the cylindrical section. Grid lines in the circumferential 
direction are generated by rotating one longitudinal plane at equally spaced angles 
around the body. Note that two circumferential planes are overlapped in order 
to facilitate the enforcement of the boundary conditions when operating in this 
direction. 

The positive orientation of the circumferential direction is chosen such that a 
right-handed system is obtained. In the present implementation of the code, this 
means clockwise for an observer looking at the body from upstream. This is very 
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Figure 41: General three dimensional view of the hemisphere-cylinder grid. 
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Figure 4.2: Details of hemisphere-cylinder typical grid planes. 
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important to ensure that anywhere on the body surface, or actually in the whole 
flowfield, the mathematical normal to any infinitesimal region is positive pointing 
outward. This becomes a crucial issue when computing particle traces with the 

graphics post-processing programs used to obtain the flow visualization figures that 
will be presented below. 


4.3 Steady State Results 

As initial tests for the code, we started with subsonic freestream cases. It is clear 
that low subsonic cases should converge faster to a steady state, because forward- 
going signals would propagate faster. However, one should be careful not to ap- 
proach the incompressible limit. Essentially the problem is that we are using a 
compressible Navier-Stokes formulation, and the system of equations becomes ill- 
behaved when the incompressible limit is approached because the energy equation 
becomes redundant. Hence, for all the subsonic cases run here, care was exercised 
to ensure that some compressibility effects were present and the above described nu- 
merical problems avoided. Subsequently, one low supersonic case was run to verify 
the capabilities of the code for supersonic problems. 

The first case run for the hemisphere-cylinder configuration was at a Mach num- 
ber of 0.5, for zero angle of attack and a Reynolds number of 1.5 million based on the 
cylindrical section diameter. Pressure coefficient contours for the converged steady 
state solution along two opposing longitudinal planes can be seen on Figure 4 3 
Mach contour, and density (nondimenaionalized with respect to freestream) con- 
tour. are al*> shown in the same figure. In this case it really does not matter whether 
these are side, or top, views of the body, because the solution is axisymmetric. 

Since experimental results could not be found for this case in the literature, the 
present computations were compared to results obtained from another finite differ- 
ence code, namely the F3D code^ 68 * 53 J generously made available to the present 
author by Ying * . Both codes have a Navier-Stokes formulation, and the lam- 
mar flow option was used. The F3D computations were performed on a Cray 2 
supercomputer, but the grid system was essentially the same as that used for the 

Ying, S.X., personal communication, September 1586. 
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(a) Pressure coefficient contours. 



u %* S « “ 

(&) Mach number contours 


-u» -*-• 


(c) Density contours. 


Figure 4.3: Flow solution for hemisphere-cylinder at M m - 0.5 and a * 0* 
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Figure 4.4: Presatire coefficient distribution on the body for hemisphere-cylinder at 
Moo = 0.5, a = 0*, and Rto 3 1.5 x 10* . 


Cyber computations in order to make the two calculations comparable. Plots of 
the pressure coefficient distributions on the body for the two computations are pre- 
sented on Figure 4.4 . As can be seen from this figure, the two computations show 
good agreement. The results from the F3D code indicate a fatter expansion over 
the hemispherical part of the body, and also predict a slightly higher magnitude 
for the negative peak Cp around the hemisphere-cylinder intersection. Both com- 
putations show this negative peak on the pressure coefficient occuring ahead of the 
hemisphere-cylinder intersection. The F3D results also seem to return faster to 
the freestream pressure value over the cylindrical section of the body. Also shown 
on Figure 4.4 is the value of the pressure coefficient at the nose stagnation point 
as predicted by isentropic relations f 37 l. These relations predict Cp * 1.06 at the 
stagnation point, and both calculations agree well with that value. 

The next hemisphere-cylinder case studied involved a fli ght Mach n umb er = 
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0.6, zero angle of attack, and the Reynolds number Re = 490000 based on the 
reference diameter. The Reynolds number in this case was chosen to match the 
experimental results presented by Hsieh^ 69,70 \ whose wind tunnel static pressure 
distributions for a hemisphere-cylinder are in the Mach number range of interest 
for this work and therefore were used for these benc hmar k computations. Com- 
putational values of pressure coefficient, Mach number, and density contours for 
this case are shown in Figure 4.5 . It can be seen from Figure 4.6 , which shows 
pressure coefficient distributions on the body for both the present computation and 
experiment^ , that the computed results agree well with the experimental ones. 
The only region along the body where the computations seem to have some dif- 
ficulty in following the experiment is on the recompression side of the expansion 
region around the hemispher-cylinder intersection. However, it is well known that 
expansion regions are particularly difficult to predict accurately. Considering that 
this is such a coarse grid system it is fair to say that the code is doing a good job 
in this case. It should also be noted that the negative peak Cp is being very well 
predicted, both in strength and location. 

Finally, it is important to point out that this a fully subsonic flow, in other 
words, there are no supersonic pockets in it. Although the ultimate objective of 
this work is to investigate transonic flows, it was felt that the grid system was too 
coarse to try to capture transonic shocks and all the possible complexity associ- 
ated with them, such as shock-boundary layer interactions and flow separation due 
to those interactions. For this reason, no transonic cases were analyzed for the 
hemisphere-cylinder configuration, and instead a low supersonic case is considered 
next. Tr ansonic flow examples were studied for the hamme rhead configurations, in 
which case much more refined grids were created. 

Computational results for a freest ream Mach number Moo = 1*5 case can be 
s een in Figure 4.7. In this example the Reynolds number was 1.386 mill ion, based 
on the reference diameter, and again zero angle of attack was considered. A plot 
of the bow shock location, calculated based on pressure gradient results, is shown 
in Figure 4.8 . It should be pointed out that in Figure 4.8 (b) we actually have 
a complete shock surface. It was plotted only as lines in order to let the body 
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MOJO. POSITION (s/ft) 


Figure 4.6: Pre>ure coefficient dietribution on hemb^herecylinder »t M» = 0.6 
and zero &ngle of >tt >c k . 
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Figure 4.7: Computed flow solution for hemisphere-cylinder at — 1.5, a = 0°, 
and Reo = 1.386 x 10 a . 
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also be seen behind the shock. The shock standoff distance can be calculated as 
—0.57 , nondimensionahzed with respect to the radius of the cylindrical section. 
This result compares well with the value of -0.60 for the shock standoff distance 
for a hemisphere-cylinder obtained from shadowgraphs, presented in Reference I 8 
It is clear, however, from the contour plots in Figure 4.7 that the computations are 
not capturing such a crisp, well-defined shock as Figure 4.8 suggests. The most 
obvious reason for this failure is the coarseness of the grid. But the main point 
we intend to demonstrate with these results, namely, that the code is capable of 
performing well for supersonic freestream conditions, is still valid. The results are 
in the correct range despite the fact we have used a coarse grid. 

Although the steady state cases studied for the hemisphere-cylinder configura- 
tion cannot be considered extremely difficult problems for many existing computa- 
tional fluid dynamics codes, the results shown provided enough confidence in the 
present code that we moved on to more complex examples. The steady state so- 
lution described above for a freestream Mach number A/* = 0.6 was used as the 
initial solution for the study of a forced unsteady case, which will be described next. 


4.4 Pitching Oscillations 

As a first actual unsteady test of the code, the hemisphere-cylinder vehicle was 
subjected to a rigid body sinusoidal pitch oscillation in a flow with = 0.6 . The 
nondimensional frequency of the oscillation was taken as 0.4, and the half-amplitude 
angle of the oscillation was $ 0 = 5° . The pitch axis was considered to be at the 
computational downstream boundary, which in this case is about 10 body diameters 
from the nose. The steady state solution previously obtained at this Mach number 
was used as initial condition for the unsteady calculation as previously mentioned. 

It should be pointed out that the frequency of the oscillation is rather large, 
perhaps two to four times larger what one should realistically expect to find in flight. 
Since expenr- eital, or computational, results could not be found in the literature 
for the type of oscillatory motion we were interested in, regardless of the value of 
frequency considered, the high value of dimensionless frequency was chosen mainly 
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Figure 4.9: Angular amplitude of centerline deflection for forced sinusoidal pitching 
oscillation. 

to allow for a fester computational turnaround time. The reader should observe 
that, for the cases considered here, the maximum allowable time step was completely 
determined by the aerodynamics. The body motion, however, is determined by the 
product UAt . It is important to note that, despite the coarse gnd being used, 
the computational costs of these unsteady calculations is not negligible. Since the 
idea behind these computations was primarily one of checking on the code, a higher 
frequency would still allow for some qualitative study of the performance of the 
code, while permitting a fester turnaround. 

The deflection of the body centerline, versus time, is shown in Fig- 

ure 4.9 . Time histories of the running normal force, i. e. , the normal force per unit 
of udai length for two axial stations along the body, can be seen in Figure 4.10 . 
It is important to point out that the “time” on both figures is the nondimensional 
time The distribution of the running normal force along the body is shown in 
Figure 4.11 for several instants of time during the oscillation. Unfortunately no 
experimental or computational results could be found to compare with the present 
calculations, but they seem to produce the kind of behavior that could be expected 
for this forced oscillation. 
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time 


Figure 4.10: Unsteady aerodynamic load on a hemisphere-cylinder undergoing a 
sinusoidal pitching oscillation with U = 0.4 (Moo = 0.6 , a = 0* , 0o = 5°) • 
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(a) Distributions for 0 = 0* (upward and downward motions) a nd for 0 = 5°. 



trior petition 

(6) Distributions for $ = 2.5 , 0. , and -2.5* , downward motion. 


Figure 4.11: Comparative plots of unsteady aerodynamic load distribution on hemi- 
sphere-cylinder in pitch oscillation = 0.6 , a = 0* , ZJ = 0.4 , S Q = 5*). 
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From Figures 4.10 and 4.11 one can see that the unsteady aerodynamic loads 
on the body lag the motion with a phase shift of almost 90°, and that the influence 
of the freest ream angle of attack at any instant along the oscillation is very small 
in determining the total load at that time instant. For example, from Figure 4.11 
we see that at 9 — 0° in the downward stroke the load is close to its maximum 
value, and at the “top” of the oscillation, 9 = 5°, the load is actually very small. In 
other words, the loads in phase with the an gul ar displacement are s ma l l compared 
to the ones induced by body oscillating velocities. Figure 4.10 also shows that the 
phase lag decreases as we move from the forebody to the aft portion of the body, 
as one might expect since the pitch axis is located at the downstream boundary. 

The reason for such a large phase lag in the aerodynamic forces, when compared 
to the body motion, is easy to understand in the light of the very high value of 
nondimensional frequency. Actually, in this case, the velocities induced by the 
motion on the forebody are comparable to freestream velocities. Since we start 
with a steady flow at time zero, we expect that the flow field will settle down 
to its steady-state, simple harmonic limiting behavior after some time. The high 
frequency also explains why these initial transients die out quickly in this case, 
which can be seen from Figure 4.12 , where the running normal force distribution 
along the body is shown for the same pitch angle position at two different cycles of 
the oscillation. The curves almost coincide, indicating that after only a half-cycle 
of oscillation the initial transients have already died out. 

In summary, for such a heavily forced oscillation, one should expect to see the 
aerodynamic loads being determined primarily by the motion itself, and the present 
computations indeed display this feature. Although the unsteady computations 
could not be 'validated against any experimental results, the correct expected be- 
havior is reproduced by the present calculations. The main objective of the analysis 
of the hemisphere-cylinder configurations can be considered as accomplished. In 
the following sections we will turn our attention to the primary focus of the present 
work, which is to apply the unsteady aerodynamic formulation to the analysis of 
aeroelastic problems in the transonic flight regime. 
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Figure 4.12: Unsteady aerodynamic l oad distribution on a pitching hemi- 

sphere-cylinder at B m 0* (downward motion) for the 1st and 2nd cycles of os- 
cillation. 



Chapter 5 

A General Hammerhead Payload 
Problem 


5.1 Preliminary Remarks 

The configuration under consideration in this section does not attempt to reproduce 
the geometry of any particular vehicle, but it has all the major characteristics one is 
likely to encounter in a hammerhead payload. For instance, it has a blunt nose and 
the payload diameter is larger than the adjacent boosting stage (which implies the 
existence of some form of boattail on the forebody). There is a flare region which is 
bound to produce a good portion of the normal force on the vehicle for *tth 1 1 angles 
of attack and that can possibly be subjected to large unsteady forces due to flow 
separation in the boattail region. 

It is clear that the study of a configuration which does not correspond to any ex- 
isting vehicle makes it difficult to compare or validate whatever results are obt aine d 
For aeroelastic analysis the problem is compounded, because not only external ge- 
ometrical similarity is necessary but also mass and stiffness similari ty would be 
required in order to compare any results correctly. However, at this point in the 
course of the research, we had not found ac^al data for an existing vehicle in the 
literature. We considered that it was worthwhile to exercise the method just to 
show what might be possible after more realistic data could be obtained. 
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5.2 The Grid System 

The computational mesh about the chosen hammerhead payload was also gener- 
ated by algebraic methods. The grid has 105 points in the longitudinal direction, 
66 points in the normal (actually nearly normal) direction, and 38 points in the cir- 
cumferential direction. Once a longitudinal plane of the grid is generated, this plane 
is rotated through 360° to create the full three dimensional grid. The rj-direction 
in this case is being called nearly normal because, in the boattail and flare regions 
of the grid, the rj- lines do not intercept the body surface at 90° angles. 

A general three dimensional view of the body and the grid is shown in Fig- 
ure 5.1 . A typical longitudinal plane of the grid can be seen in Figure 5.2(a) , 
and details of the forebody and flare regions are shown in Figures 5.2(6) and (c) . 
The latter illustrates the rj-lines not intercepting the body surface at right angles 
in the boattail and flare regions. This is not the ideal way of creating a computa- 
tional mesh, since in principle we want the normal coordinate lines indeed coming 
in normal to the body. A small numerical error is expected due to the fact that 
this is not precisely true. However, numerical computations of boattail flowfields by 
Deiwert^ 57 ! indicate that this does not cause enough numerical problems to justify 
a more elaborate grid generation scheme, if the boattail (or flare) does not have 
a very steep slope. Moreover, when performing aeroelastic analyses, the grid will 
deform as the computation proceeds, and it can be very difficult to ensure that r? 
grid lines do not cross over each other if the grid generation scheme gets to be too 

complex. 

Since true transonic solutions are attempted, it is important in this case to have 
a better clustering of the grid points along the body (in the longitudinal direction) 
in regions where high gradients are likely to occur. For this reason, one parameter 
hyperbolic tangent grid stretching t 561 is used to cluster the grid points around 
the hemisphere-forebody cylinder and the flare-afterbody cylinder intersections. In 
order to improve the resolution on the upstream centerline region, one-parameter 
hyperbolic tangent grid stretching is also used to cluster grid points around the 
upstr eam centerline. To avoid wasting grid points on the downstream part of the 
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Figure 5.1. Three dimensional view of general hammerhead payload configuration 
grid system. 
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(a) Complete plane. 
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( b ) Detail of forebody region. 
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(c) Detail of flare region. 


Figure 5.2: Typical longitudinal grid plane for hammerhead payload. 
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afterbody cylinder, where not much is happening, a 8.5% exponential stretching is 
used to increase the grid spacing towards the downstream boundary. In the nor mal 
or nearly nor ma l, direction a 16% exponential gnd stretching is used everywhere in 
order to cluster grid points near the body for capturing viscous effects. 


5.3 Initial Aerodynamic Solution 

As previously explained, before the structural-dynamic equations are turned on, 
an initial aerodynamic solution must be obtained at the desired flight condition 
under the assumption that the body is rigid. This is done by initializing all the flow 
variables with freestream conditions, imposing the appropriate boundary conditions 
for the problem, and then letting the aerodynamic solver march the solution in time 
until some steady state condition is achieved. An important caution, however, is 
that one should be aware that, for these transonic freestream conditions, a true 
steady state condition is not necessarily obtained. It is possible that the aerodynamic 
solution by itself is unsteady. In these cases what is being searched for in this initial 
phase is is a correct possible solution for the flowfield around the body at that 
particular flight condition. 

For the examples considered here, however, true steady state solutions were 
obtained in this initial phase. The steady state flow solution was calculated over 
this geometry for a freestream Mach number Af* = 0.85, an angle of attack a = 6°, 
and the Reynolds number Re = 1.26 x 10 e (based on the diameter). The flow, or 
rather the boundary layer, was considered turbulent, i. e. , the turbulence model 
was turned on. Pressure coefficient contours, Mach number contours density 
contours for leeside and windside are shown on Figure 5.3, where essentially we are 
looking at a side view of the body. The expected flowfield features are reproduced 
by the computation. For instance, the expansion regions around the hemisphere- 
forebody cylinder intersection and around the flare-afterbody cy lin der intersection 
are clearly shown on the figures. A mild expansi on around the forebody cylinder- 
boattail intersection can be seen, and the compression region on the face of the 
flare, mainly on the windside, is also very well defined. 


63 



rw rt- 


ORIGINAL PAGE IS 
OF POOR QUALITY 

CHAPTER 5. A GENERAL HAMMERHEAD PAYLOAD PROBLEM 



(b) Mach number contours. (c) Density contours. 

Figure 5.3: Flow solution about a hammerhead geometry at Mcc — 0*85 , constant 
a = 6* and Reo = 1.26 x 10* (side view). 
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All these features are even more apparent from Figure 5.4 , which shows leeside 
and windside pressure coefficient distributions on the body. As mentioned, these 
pressure distributions seem to have all the characteristics that would be expected, 
except that a more clearly defined shock was expected to be seen on the leeside of 
the forebody. Although we do not have experimental data with which to compare 
these pressure coefficient distributions, it is worth noting that xfhe magnitudes of 
the negative peak Cp on the hemisphere-forebody cylinder intersection are in the 
correct range for this flight condition. 

An interesting point that was discussed previously concerning the downstream 
boundary conditions can be clearly understood from Figure 5.4 . One can see that, 
on the afterbody cylinder, the pressure coefficient has a constant value different 
from zero for a good portion of this body section; it then smoothly drops to zero 
at the downstream boundary. Of course, there is no surprise that Cp is zero at 
the downstream boundary, because we are enforcing that as boundary condition. 
However, the existence of a fiat region in the Cp distribution ahead of the down- 
stream boundary evidences that there is an error in the boundary condition. If the 
computational domain is extended further downstream, perhaps with the inclusion 
of a body base, the flow would converge to a nonzero value of pressure coefficient 
at the body location where the boundary condition is currently being enforced. On 
the other hand, the existence of this flat region is also an assurance that w ha tever 
errors are being produced at the downstream boundary are not propagating far 
enough upstream in order to influence the region of real interest in this case, which 
is mainly the forebody. 

It should be clear that the amount of data generated in these three dimensio nal 
computational solutions is very large. The study of scalar flow variable contours 
may prove inadequate to understand what is really happening. To this end, particle 
traces are very useful. In particular, particle traces restricted to the first computa- 
tional rj - plane away from the body amount to computer generated otUflow pictures 
and are almost indispensable to understand the topology of the flow. It should be 
pointed out that the latter cam also be interpreted as plots of the skin-friction lines 
on the body. We shall not go into any detailed discussion of these flow topologies, 
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(6) Windside. 

Figure 5.4: Pressure coefficient distribution along the body for a hammerhead pay- 
load at Moo = 0.85 , a = 6* , and Reo — 1.26 x 10* . 
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except to point out lines of separation and/or reattachment which will help us un- 
derstand the solutions that are being obtained. The reader interested in detailed 
discussions of topological flow structures is refered, for example, to the works of 
Dallmann^ 71, Kaynak, Holst and Cantwell ^ 73 ^, and DeiwerJ 74 ^ . 

Figure 5.5 shows these computer generated oil-flow pictures for this configuration 
at the converged steady state solution. One can clearly identify a separation region 
on the cylindrical forebody right after the hemisphere-cylinder intersection, and a 
node of separation can be seen on the lee generator just behind this intersection. 
There is another separation line on the boattail, which indicates that even the 
windside experiences some flow separation in that region. The line of reattachment 
aft of the boattail is also clearly defined in the figures. All those cases can be 
considered mild separations, in the sense that the regions of reversed flow are rather 
limited. Figures 5.6(a) and 5.6(6) show velocity profiles for the leeside in the regions 
of separation on the forebody cylinder and boattail, respectively. The reversed flow 
directions are apparent, indicating the backflow condition in the separation regions. 

The convergence to steady state was rather slow, which is expected for a tran- 
sonic flow condition. To achieve convergence, 4000 time steps (iterations) were 
needed. It should be mentioned that, for the 105 x 66 x 38 grid being vised, each 
iteration takes approximately 10 CPU seconds in a CDC Cyber 205 computer. The 
overall system time, however, is a little higher than what would be estimated from 
this figure, because the database is not core contained. Data must be shifted back 
and forth from disk. As explained before, the database is structured in a pencil 
format, and the metrics of the transformation are what is kept in outside disk files. 
When operating in a particular direction, the metrics for that direction must be 
read in from those files. Due to the use of asynchronous I/O, this can be done with- 
out degrading very much the performance of the code for steady state problems. In 
the present example, the CPU time represented approximately 60% of the overall 
system time. 


67 



CHAPTER 5. A GENERAL HAMMERHEAD PAYLOAD PROBLEM 



(a) Side view. 



(6) Top view. 


Figure 5.5: Computer generated oil-flow lines for general hammerhead payload at 
Moo — 0.85 and a = 6° . 
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(a) Forebody cylinder region. 



(b) Boat tail region. 


Figure 5.6: Velocity profiler on the leeeide in the region* a £ flow reparation (tide 
view). 
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5.4 Aeroelastic Analysis 

The steady state aerodynamic solution described in Section 5.3 will be used as 
the starting flow for the aeroelastic cases studied here. Since this solution was 
calculated for a rigid vehicle, the airloads are not the correct actual loads at a 
deformed equilibrium position for the elastic vehicle. This fact provides a way of 
introducing the initial perturbation to start the oscillation, which will be adopted 
in the present work for all cases where the freestream angle of attack is different 

from zero. 

The structural information necessary for the present analysis consists of some 
normal mode shapes and their corresponding natural frequencies, the structural 
damping coefficients associated with each mode, and the vehicle mass distribution or 
the generalized masses associated with each mode. As previously discussed, we have 
no actual data for an existing vehicle that matches this configuration. Therefore, 
the approach foUowed consisted of estimating the necessary properties from what 
could be gathered from the literature. For instance, three structural modes were 
employed in this case, and the chosen mode shapes resemble those presented by 
Woods and Ericsson ™ . 

The theoretical analysis of the numerical stability of the scheme selected for the 
solution of the aeroelastic equations showed that the scheme is numerically nondis- 
sipative. This is a very important issue, because we want to make sure that the 
numerical scheme is not stabilizing physically unstable aeroelastic solutions due to 
numerical dissipation introduced by the method. The first example here undertakes 
to address this problem. Essentially, we want to demonstrate numerically what the 
Unear stabiUty analysis already predicted, i. e. , that the numerical method is not 
intr oducing any numerical dissipation (or instabiUty) and that a pure sinusoidal re- 
sponse can be captured. In this first case run, therefore, all the structural damping 
coefficients were set to zero, and a very small value of dynamic pressure was con- 
sidered. The s mall dynamic pressure is chosen so as to minim ize the damping effect 
of the airstream. The response, in terms of the generalized modal displacements in 
the second and third modes, is shown in Figure 5.7 . The results indeed confirm 
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Figure 5.7: Response for zero structural damping and very low flight dynamic 
pressure. 
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what was expected, since an undamped response is obtained. Although this was the 
anticipated result, it is important to see that all the complex algorithms developed 

are able to reproduce such a simple phenomenon. 

The other cases analyzed considered more realistic values of dynamic pressure 
and also some nonzero values of the structural damping coefficients. The proce- 
dure followed was to keep the flight Mach number and the angle of attack con- 
stants, and vary the dynamic pressure. The structural parameters were, of course, 
kept constant. For instance the values of structural damping coefficient used were 
Ci = 0.0010 , Ci = 0.0018 and Cs = 0.0036 , for the first, second and third modes, re- 
spectively. A typical vehicle response for an intermediate value of dynamic pressure 
can be seen in Figure 5.8 , which shows the response in each of the three gener- 
alized coordinates. By comparing the magnitudes of the responses for the three 
modes, it is clear that the overall vehicle response is dominated by the first mode 
displacements. Figure 5.9 makes this point even more clearly by showing the total 
deflection at the nose of the vehicle, which can be seen to be composed of the first 

mode deflection plus a small higher frequency influence. 

The usual way in which the influence of the aiistream is determined consists of 
c al c u lating the damping coefficient for the run with the air on, in other words, 
for dynamic pressure different from zero, and comparing the rate of decay of the 
motion with the pure structural damping coefficient. If the former is smaller than 
the latter, we have a situation where the presence of the flow is destabilizing. Of 
course, in order to obtain an unstable condition this influence must be large enough 
to overpower the structural damping and cause the amplification of initial pertur- 
bations. As one can see from Figure 5.8, the response in all three modes is damped 
in this case. The response on the first mode is only slightly damped, the second 
mode shows a little faster decay, and the third mode decays the fastest. A closer 
analysis of the damping coefficient reveals that, in this example, the first and sec- 
ond modes are rendered even more stable by the airstream. On the other hand, the 
damping coefficient for the third mode, at the present dynamic pressure level, is a 
bit smaller than the structural damping Cs = 0.0036 . One therefore cncludes that 
the freestream is feeding energy into that mode’s oscillation. 
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Despite the fact that them « time donmin amtly*., the remit, are beet sum- 
mlri2M bjr a root locus plot, which is shown in Figum 5.10 . In this plot the arrows 

indicate the direction of increming dynmnic pressure, which » the 

"he analysis. The abecssa is the re^ part of the «roefastic root whrch ts a 

measure of the rate of the decay of the collation in ««h mode, formed * 
net of the damping coefficient at that particular dynamic premure tune, the natuml 
CC for that mode. The ordinate is the una P nmy par, of the aeroelasticroo. 
which is the frequency of the response of that mode at the dynanuc P»"™ ’ 
errf Although one should note that the malm on the two axe, are very drfferent^ 
it can be seen from the plot that the frequencies reuuun apprmnmately unc g 
throughout the whole range of dynamic pressures considered. 

All the eases analysed for this hmnmerhmd shape were aeroefastrcally stable, 
although all the mode, showed an initial tendency of going toward, the unstable 
side for very smJl value, of dynamic pressure. This tendency wm ^ 
as the dynamic prewrure was increased further, such that it was 
soy of this behavior from the time plot for the tot mid second modes. For the 
third mode, however, the reversal of damping at low freestream dynamic pressure 
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Figure 5.10: Aeroelaatic root loci for general hammerhead configuration at 

00 0.85 and ot — 6 . The parameter varied is flight dynamic pressure* 
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is clearly visible from Figure 5.10 . After this initial trend was passed, what we 
wiU call a pare damping behavior, for the lack of a better description, was observed. 

It is characterized by a simple increase of the damping in each mode as the flight 
dynamic pressure is raised. The frequencies remain approximately unchanged. This 
is not typical aeroelastic behavior for conventional flight vehicle configurations with 
wings and tails. However, it should be noted that the natural frequencies considered 
for this analysis are quite high, which means that we are probably assigning stiffness 
values that are higher than they should be and explains why no flutter is observed. 
Moreover, as already mentioned, this configuration does not correspond to any 
existing vehicle. It is, therefore, very difficult to try to correlate these results with 

some expected behavior. 

Finally, the issue of the computational costs of these aeroelastic solutions should 
be considered. For the 105 x 66 x 38 grid which was used here, each aeroelastic iter- 
ation takes about 12 CPU seconds in a CDC Cyber 205. This represents an increase 
of 20% over the CPU time per iteration for steady state aerodynamic calculations. 
Almost all the additional time is spent in regenerating the computational grid once 
a new body deflected position is determined from the solution of the structural- 

dynamic equations. 

One does not necessarily have to regenerate the complete grid at every time step, 
but some form of reshaping is required in order to account for the deformation of 
the body. Some authors, see for instance Steger and Bailey 1 11 , prefer to use some 
form of shearing transformations to account for the grid deformation, and avoid 
solving the grid generation equations at every time step. In the present approach, 
since an algebraic grid generation scheme is bang used, we regenerate the whole 
grid. The grid points on the body surface keep their relative position with respect to 
the body centerline at every axial station, and the far field boundary is kept fixed. 
The remainder of the grid, which is all the interior part of the grid, is interpolated 
between those ends using the same algorithm that initially created the whole grid. 
Note that the structural-dynamic equations are solving for the deformation of the 
body centerline, which means that most of the grid displacement is going to occur 
close to the body. This is consistent with what is physically happening, in the sense 
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that the body motion should cause perturbations close to itself and the far field 
should remain undisturbed. 

The amount of time consumed by the solution of the stnictural-dy namic equa- 
tions IS truly negligible when compared to the overall CPU time. The point that 
should be considered, though, is that modal analyris is being used. As previously 
mentioned, this is a very powerful technique since a whole range of effects can be 
considered. However, one may want to replace the modal superposition approach 
y, or example, a finite element representation of the vehicle. In such a case the 

amount of time required for the stnictural-dynamic solution is bound to increase 
considerably. 

The overall system time per iteration for aeroelastic analysis, however, suffers 
much more than the CPU time when compared to a pure steady state aerodynamic 
solution. As discussed in the previous section, the major difficulty arises from 
t e fact that the datable is not core contained. In this case, since the grid is 
recalculated at every iteration, so are the metrics of the trenrfbrmation. Essentially 
we me doubling the amount of I/O per iteration, since now w. have to output 
« metrics re they are calculated mid then read them back as we operate in each 
dmection. Although asynchronous I/O is still bring used, the performance of the 

d of CPU time to overall system time is only 

«ound 30% Tbs means that, TO % of the time the job is running, it is simply doing 

mpu /output Of course, there kinds of statistics are very mrehin. dependent, mu! 

Z m e T >h ‘ CD ° CybW 205 ' “"*• « —M >*e to avoid all 

Y ^ v *“* ° f Pr ° bl ' m ” “ dKjin * Wre this could certainly be 

accomplished by uring a computer such re the Cray 2 . Nevertheless, all the effort 

spent m maloug the algorithm suitable even for “smaller” computers is important 

becaure show, that the method could be used even by there without access to a 
system of the size of the Cray 2. 
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Chapter 6 

Analysis of an Atlas- Able IV 
Configuration 


6.1 Initial Considerations 

As previously mentioned, in the early 1960’s some aeroelastic problems were ob- 
served on launch vehicles with hammerhead payload configurations when passing 
through the transonic regime. One of those configurations was the Atlas-Able IV, 
in which case the problems observed in flight were later traced back to something 
resembling flutter. The mechanism driving the instability was associated with the 
lags in the aerodynamic forces caused by phenomena induced by the hammerhead, 
which would cause the airstream to do positive work on the oscillating vehicle. 

The ideal case to study here would be one where both steady aerodynamic data 
and detailed aeroelastic results were available. Unfortunately, this is certainly a 
difficult combination to find in the literature. At the time, the Atlas-Able IV con- 
figuration represented the best test case we could find, despite the fact that aeroe- 
lastic information is somewhat limited. Geometrical data for this configuration are 
available, at least in nondimensional form* , and steady aerodynamic pressure 
distributions from wind tunnel tests are also published* 7 1 . The structural data, 
however, are rather sketchy because part of this information is still classified 

* Ericsson, L.E., personal communication, April 1987. 
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For instance, exact values for frequencies, structural damping coefficients, and ve- 
hicle mass distribution were not accessible from the same sources that provided 
geometrical and steady aerodynamic data. 

In the present work, the approximate range of frequencies to be considered, the 
approximate form of mode shapes, and an estimate of reasonable structural damping 
coefficients were obtained from Woods and Ericsson and Ericsson * . The mass 

distribution was estimated from data available for Atlas boosters and provided by 
Gen. Dynamics. Steady aerodynamic wind t unne l data, which was used to validate 
the initial aerodynamic solution, was available from Graham and Butler^ 75 ! . 

Other steady and unsteady aerodynamic data were found in the literature for 
similar hammerhead configurations. For instance, Coef 78 ^ presents the steady and 
unsteady aerodynamic pressures on an Able V payload model, and Robinson et 
af-^ 7 ^ study the dy nami c response of some hammerhead models to unsteady aerody- 
namic loading. However, the Atlas- Able IV data still constitute the best information 
and allow for a better comparison of results. Finally, it should be mentioned that 
a very complete set of structural and aeroelastic parameters on two Titan/Centaur 
models was made available to the author by Henning • Unfortunately, because 
of time limitations and the high computational costs of these aeroelastic solutions, 
it was not possible to perform aeroelastic simulations for thes e cases. 

As we have done in the previous cases studied in this work, the computational 
gnd for the Atlas-Able IV configuration was generated using algebraic methods. A 
three dimensional view of the body and grid in this case can be seen in Figure 6.1 . 
The grid has 105 points in the longitudinal direction, 66 points in the normal direc- 
tion, and 38 points in the circumferential direction which A gain is a periodic type 
mesh going 360* around the body. In this case, the body is composed of a ellip- 
soidal nose, a cylindrical forebody section, a boattail, and the cylindrical afterbody 
section. .-The ratio between minor and major axes of the ellipse is approximately 
0.605 . The reference length adopted in this case is the diameter of the cylindri- 
cal afterbody section, and the diameter of the cylindrical forebody Taction is 1.5 
reference lengths. 

A typical complete longitudinal plane of the grid is shown in Figure 6.2(a) , 

Ericsson. L.E., personal communication, April 1987. 

+ Henning, T., personal communication, May 1987. 
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Figure 6.1: Three dimensional view of Atlas-Able IV configuration grid system. 
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details of the nose and the boattail regions can be seen in Figures 6 . 2 ( 6 ) and (c) , 
respectively. In the longitudinal direction, 29 points are used in the nose region, 
12 in the cylindrical forebody region, 27 in the boattail, and 37 in the cylindrical 
afterbody region. One parameter hyperbolic grid stretching^ ^ techniques are used 
in the ellipsoidal nose region to cluster grid points towards the upstream centerline 
and towards the ellipsoid-cylinder intersection. An equally spaced grid is used over 
both the forebody cylinder and boattail sections, and a 6.98% exponential grid 
stretching is used in the cylindrical forebody section in order to gradually coarsen 
the gnd as we move towards the downstream boundary. 

The normal direction (or 17 -direction) would again be more properly called a 
nearly normal direction, since r?-lines do not intersect the body surface at exactly 
right angles on the nose and boattail regions. Since these lines meet the body at 
angles that are very close to 90° , we will accept the small error being introduced 
by this simplification, as already discussed in the previous chapter. To ensure 
proper capturing of viscous effects in the normal direction, a 16% exponential grid 

stretching was used everywhere in this direction in order to cluster grid points close 
to the body. 


6.2 Configuration at Angle of Attack: Steady State 
Results 

The initial aerodynamic solution involving the Atlas-Able IV was calculated for 
a = 0.85 and a = 6 ° flight condition. The Reynolds number considered was 
Re = 1.2637 x 10* , based in the reference diameter. It should be mentioned that this 
number was chosen to match Graham and Butler’s experimental conditions [ 751 . 
The boundary layer was considered turbulent for most of the calculations performed 
here. Starting from freest ream everywhere and imposing the appropriate boundary 
conditions, we allow the flow to evolve to a converged solution. 

A side view of the pressure coefficient, Mach number, and density contour plots 
around the body for the converged solution can be seen in Figure 6.3 . These plots 
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(a) Complete plane. 



■ 


(6) Detail of nose region. 



Figure 6.2: Typical longitudinal grid plane for Atlas-Able IV configuration. 
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(6) Mach number contours. (c) Density contours. 


Figure 6.3: Flow solution about an Atlas-Able IV configuration at Moo — 0.85 and 
Qt = 6° (side view). 
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xn 



Figure 6.4: Pressure coefficient distribution on the leeside of an Atlas-Able IV 
payload (Af*, = 0.85 , a = 6°). 


give a good idea of the overall flowfield appearance, despite the fact that we are only 
seeing leeward and windward planes. The Mach number contours seem to indicate 
a sizable separated region on the leeside, which is evidenced by the bubble sort of 
behavior of the contours over the forebody cylinder and boattail regions, where the 
local velocity magnitude (or Mach number) increases, then decreases back to zero, 
and finally increases again up to the freest ream condition as we move outwards from 
the body towards the farfield. This is very typical of reversed flow regions, and it 
will be investigated in more detail later when discussing particle trace results. 

Taking advantage of the fact that for this configuration there are experimental 
results available, Figure 6.4 shows a comparison of the computed pressure coefficient 
distribution on the leeside of the body with wind-tunnel measurements by Graham 
and Butler 1 .It can be seen that the calculated Cp distribution follows the trend 
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of the experimental results well, but the strength of the shock is not bring accurately 
captured by the computation. The location of the shock seems to be correct, but 
apparently the flow doe, not expand as much as , h , experiment indicates it should. 
The result u a weaker shock. It is clear that the computation* shock is spread over 
a few gnd points, which is typical of centrally differenced finite difference scheme, 
The computat.ons indicate that the flow separate, right after the shock, which 

f "IT'”" 11 Wi “‘ ScUimn ph0 ‘ 0 *»P>“ liable from Reference' 75 ' 
1 fllght COn f‘° n ' AS f ° r * h « pre ““« coefficient distributions, the sepa- 
rated regmn » evidenced by the somewhat flat Cp distribution. The success of the 

present computat.on of numerical value, of Cp over the separated region is also 

.. h f P '1" 1 ' h ° Ugh C ° mPUtati °“ “ pture lhe comct trends in the Cp 
“ kITu mUS ‘ h* POinted thAt Terence in the length of 

“ "* , ° n bel ’ reen the *«*“*ry provided by Eriowm * and the one 

c P Z^«t "'T ^ BU ‘ ler[75 ' ' ^ ll “ CT *“ * *«•« fore body 

" j. and ’ con5w l u «*ly. a slightly longer boattail action such that the 

ovmall payload length is the same. This make, the Cp comparison, over the boattail 

bTffifLI T ^ dUe l ° ' he Pr0b,OT> "*** ““ strength it would 
, any event, to compare result, downstream of the shock The prob- 

: b :irrr by * he ** ,hit * th - ~ — •*>. 

tobnlenc. modelhng becomes even more important, and this may be afljfog the 

A compmwon of the premure coefficient distribution, on the body for other Ion- 
Cw^d^r b ;“” " 6 5 ' ^ Cp — over the 

mTrh ,ortb ~ otb — • ~ zszzz 

over the 11^1, T ’™* 1 »»“*. »•“ «e slightly undemxpanded 

ofamrJl^T^ oo^forebody cylinder intemection. A new feature that can be 

a second d” ~ MU)t ’ “ ^ «“» “ ‘ -iceable tendency 

The won for°ilt' C ' P U ° Und th * f ° r ' b ° dy c > ,| “drr-boattaU intersection, 
he won for that » qurte clem, since it i, wonable to expect the sepamtion to 

be less severe on the aide of the vehicle, or on its windside, tZT it is on“ I 
Ericsson, L.E., personal communication, April 1937. 


85 




86 


ORIGrii'fAL PAGE fS 

OE ft)OR QUALITY 




CHAPTER 6. ANALYSIS OF AN ATLAS-ABLE TV CONFIGURATION 


This second dip is merely • result of the flow reaction due to the expansion comer 
over the forebody cylinder-boattail intersection. The experimental results for the 
wmdwmd^ne exhibit what seems to be some data matter, as on. can see from 
Figure 6 . 5 ( 6 ) . Reference' 75 ' offer, no explanation for its existence. Our calcu- 
lations do not repmduc. any of this behavior, and the calculated Cp distribution 

oUows the experimental curve that would be obtained if the scattered points were 
neglected. 

Finally, it is important to stress that the compariwn. were nmde between the 
computational result, and the experimental one, at the same nominal tunnel Mach 
number. In other words, there was no attempt to try to correct for the tunnel 

from^f ’ ! th ' ' ff “* iVe M “ h nUmber “ th ' ta— be different 

e nominal one. Usually one tnes to match the vehicle lift coefficient instead 

o matching the nominal tunnel parameters. This was not done in this case, and it 

Z. n' T ° f <UsCr ' P “ d “ m the Cp result,. Them were 

me numerical difficulties, mrnnly associated with artiflcUd dissipation pammeters 

pZ I t ° f ““ <*■— - ‘be results We would 

the next Ziom 0 ^ ^ * U “ ““ -blems until 

A good pictorial description of the flow topology cm, be men from Figure 6 6 

at tLffiZ S ? nd l ° P **” 0a ' fl0W Un " for lh * Atlm-Able IV configuration 

*“ - ^ fl ° W ■*“*- *bn«* *« ‘bo ellipsoidal 

th. wav d , C Z ° n “ “* ‘*“ ide - T1 “ ° 7 separation extend, all 

of thTfolZv WW ,h ' flOW enmewhere just downstremn 

of ‘befombody cyhnder-boattail intention. In order to help the understanding 

restricting^’ I r “"““““S ‘bn*’ «-ce ofl-flow line, am generated by 

^j Cl ' “ l ° ,h * 3<con<i Pbm' »bove ‘be body surface, lines of sep. 

oZThZ 7 f ^ wh “ »-«* •»- converge to it. On the 

her hand, lures of reattmhment me shown by empty spaces, whem all trace, go 

««y om it. Th' f gum also shows that the boatt.il is completely immemed in a 

^m of reversed flow, mrd that the flow reattaches M mewhem downs, mmn of the 
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(a) Side view. 



(6) Top view. 


Figure 6.6: Oil-flow tine, for Atla.-Able IV configuration at = 0.8S and a 
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Unrestricted particle traces are shown in Figure 6.7 for some different viewpoints, 
which help visualize and understand the overall flow topology in this example. 
Figure 6.8 shows an expanded view of the particles taking off from the body surface 
around the lee generator region. It is clear from these figures that in this case 
we have a very extensive region of reversed flow, and of course flow separation. 
Figure 6.9 complements this picture by showing velocity profiles on the lee and 
windsides over the regions of reversed flow on the body, where we can see how far 
the reversed velocities penetrate into the flowfield. This information reassures us 
of our interpretation of the Mach number contours presented in Figure 6.3 , which 
have been previously discussed. It also makes it evident that even on the windside 
there is a region of separated flow for this flight condition. Finally, Figure 6.10 
shows a front view of the vehicle, with particle traces indicating that the solution 
is indeed symmetric with respect to the pitch plane in this case, which confirms 
statements previously made. 

It is very important to realize the difference between the flow solution obtained 
for the present configuration and the one for the configuration studied in Chapter 
5 . Although the present configuration is a somewhat more slender one, the extent 
of the flow separation region in the present case is much larger than on the previ- 
ous case. It is the author’s belief that such behavior is mainly associated with the 
length of the forebody cylinder region. This region is short in the present configu- 
ration, which allows for the merging of the distinct separation regions observed in 
the configuration of the previous chapter. In other words, the separation caused 
by the shock, which usually is located right after the nose region at least on the 
leeside, merges with the flow separation due to the adverse pressure gradients on the 
boat tail. This creates a very extensive separated region, which explains the kin d of 
flow topology observed in this case. Moreover, the shorter forebody also produces 
more severe adverse gradients, which again contribute to extend the separated re- 
gion. The existence of such extensive flow separation raises questions with respect 
to the aeroelastic stability of the configuration, besides being a source of conce 
with respect to buffeting loads. 

Since the grid size used in the present investigation is the same as in the case 
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Figure 6.7: Particle traces showing flow separation on the Atlas-Able IV configura- 
tion at angle of attack. 
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Figure 6.8: Expanded view of flow separation close to the leeside. 
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studied in the previous chapter, i, «. , a 105 x 66 x 38 point grid, the computational 
statistics are essentially the same as on the previous case at least as / 

state calculations are concerned. For a CDC Cyber 205 computer, ~" lon 
takes approximately 10 CPU seconds, mid somewhere around 35 “ * ™ 

step, (iterations) are required to achieve a converged solutton. The CPU 
represents approximately 60% of the overall system time, and the remainder 
time is spent in I/O because the datable is not core contamed. 

6.3 Some Computational Difficulties 

This section discus*. some of the difficulties encountered in ‘tecomputation ofthe 
aerodynamic Bow, studied here. It should be noted that we wdl he ridkmg about 
steady state calculation, in the present section, since meet of the vahdatton of 
p^y aerodynamic computational method was done for sternly state cases. It . .ho 
LC- to realise that, s.nce the use of superp-ition in a modal form has- 
certainly proved for aeroelastic analysis, the major cone™ in the pmsent work .s 
Show that the computational technique, propewed here are capableof *ProduJ 
the cormct aerodynamic phenomena nec«sary for the aeroelast.c analysts. For tbs 
reason it is very important to study these initial aerodynamic solutions, or steady 
state solutions, to ensure that the physical flow phenomena are being adequately 

C&Pt nL seen in the previous action that the cod. is doing a less £-*«<£ 
job in terms of capturing the shock strength. On. of the pcesibbexplanatto^o 
this kind of difficulty is a-ociated with the amount of artffim J <h~.pet.on bemg 
used in the calculation. A. discussed before, for central differeneeschem«i,s 
necessary to introduce some form of numerical disripatkm m order to control 
nonlinear instabilities amocieted with the aliasing back into the lower frequ«.cy 
range of the high frequency phenomena that is not supported by the mesh. The 
question is, then, how much artificial dissipation should be added. 

The approach followed in this work consisted of addin* the mnumum amomri 
of artificial dissipation that would still ensure numerical stabdrty of the solu 
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process. This minimum amount was determined by numerical experiments. When 
one starts from freestream and tries to march the solution in time up to a steady 
state condition, the time steps that have to be taken in the beginning of this process 
are very small. One way to speed up the convergence to steady state is to take larger 
time steps and simultaneously use unrealistically high values of artificial dissipation 
in order to keep the solution numerically stable. As the flow solution approaches 
the correct steady state condition, the numerical dissipation can be reduced to more 
appropriate levels where it should not interfere with the physics of the computation. 

The problem observed in the present work was that the amounts of artificial 
dissipation necessary to keep the numerical stability of the solution process were 
almost one order of magnitude higher than what is recognized in the literature * 
as the appropriate level for these computations. One possible consequence of this 
use of excessive amounts of artificial dissipation is that the shock might have been 
dissipated by them. This can be another factor explaining the features observed in 
our computational results with regard to the transonic shock capturing, where the 
flow does not expand as much as the experimental results indicate it should, and 
consequently the shock strength is not correct. 

It should be mentioned that the artificial dissipation scheme imp lem ented in 
the present code is what can be called a constant coefficient artificial dissipation 
algorithm. This means that a constant value of artificial dissipation coefficient is 
used throughout the computational domain. More recent flow solver algorithms, 
still using central differences, have more elaborate artificial dissipation algorithms, 
usually weighting the artificial dissipation coefficient with the local norm of the 
residue. The latter seems to produce better results, but unfortunately we did not 
have the opportunity to implement it in the calculations performed here. 

Another subject that presented considerable concern in the present work was the 
issue of turbulence modelling. It is understandable that this should be true, since 
the models are of an empirical nature. As we have mentioned before, the two layer 
Baldwin and Lomax algebraic eddy viscosity model was used for the present 
computations. This model was originally derived for attached or mildly separated 
boundary layers. The results in the previous section indicate that a rather massive 

Pulliam, T.H., personal communication, October 1987. 
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flow separation condition is being observed in this case. 

Actually, the problem becomes more complex by the feet that some small vari- 
ations in the model computational parameters can cause dramatic changes in the 
flow topology. The parameter varied in the present work was the distance, in terms 
of computational grid points, from the body wall that we search for the maximum 
velocity in the profile, which is an important parameter in the present model imple- 
mentation. In the results presented in the previous section, this search is done up 
to the 25th grid point in the normal direction. This was the value of this parameter 
(called edgek in the present work) that seems to produce most realistic results in 
terms of the flow topology. For values of edgek smaller than 25, the separation 
region is even larger, and for values larger than 25 the flow shows a tendency to 
remain attached over larger portions of the body. 

Figure 6.11 shows oil-flow lines for a solution obtained using edgek equal to 20 , 
i. e. , the ya r^ h for the max imum velocity in the boundary layer profile is done up to 
the 20th grid point. The freestream parameters are the same as before, Moo = 0.85 
and a = 6°. In this case there is a well defined foeus^ 1 ^ on the side of the 
body. Even ahead of the ellipsoid-cylinder intersection on the lee generator there 
is a saddle point of separation, whereas on the solution presented in the previous 
section a nodal point was observed in the corresponding position. The release of 
particles around the focal point, as identified from Figure 6.11 , produces the particle 
traces plot shown in Figure 6.12 . This figure shows how the particles are caught 
in the reversed flow region and convected upstream before reaching the forward 
flowing stream region. Details of the vortex t&king off from the body surface can be 
seen in Figure 6.13 , for the same particle trace plot shown in the previous figure. 

It is clear from these figures that the topological structures observed in this 
are very different from those observed on the results presented in the previous 
section. Moreover, the region of separated flow in the present case is larger than 
in the previous one. Although the flow reattaches just downstream of the boattail, 
as before, the separation occurs further upstream when compared to the previous 
case. In this case not only the boattail, but also most of the forebody cylindrical 
section are immersed in a reversed flow condition. Finally, it should be mentioned 
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(a) Side view. 



( b ) Top view. 


Figure 6.11: Oil-flow lines for computation with search up to 20th grid point. 
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Figure 6.12: Side view of traces for particles released around the focal point on the 
side of the body. 

that, even with such dramatic flow separation, the solution is still symmetric with 
respect to the pitch plane, as before. 

We turn our attention to the cases when the parameter is varied in the other 
direction. Figure 6.14 shows a side view of the oil-flow lines for the solution obtained 
when the search for the m^yimum velocity in the boundary layer profile is performed 
up to the 36th grid point. The flow on the forebody cylindrical section is fully 
attached in this case, and the separation region is limited to the upper portion of 
the boattail region. In terms of flow separation structure, the swirling of the flow 
around the focal point on the boattail is evident from the figure. It is clear from 
this figure that the flow on the windside never really separates, which is also in 
contrast with the results previously shown for this configuration. 

The parameter can be further increased, and when the search is performed up to 
the 45th grid point the solution is fully attached, as evidenced by the oil-flow lines 
shown on Figure 6.15 , where we are looking at a side view of the vehicle as before. 
At this point, further increases of the parameter will not change the topology of the 
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Figure 6.14: Side view of oil-flow lines for computation with search up to 35th grid 
point. 



Figure 6.15: Side view of oil-flow lines for computation with search up to 45th grid 
point. 
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flow. However, when the search for the maximum velocity on the boundary layer 
velocity profile is performed up to the 55th grid point, some numerical stability 
problems start to be developed on the nose region. In this last case, i. e. , search 
up to 55th gnd point, although a solution that looks very much like the one shown 
on Figure 6.15 can be obtained, the maximum residue could not be dropped below 
10 -3 and it was actually increasing at the point we stopped the computation. 

Perhaps the most interesting point of all this discussion, at least as far as aeroe- 
lastic applications are concerned, is that despite the dramatic changes in flow topol- 
ogy the pressure coefficient distributions on the body are not much affected by all 
these variations on the parameters of the model. Of course, there are some alter- 
ations in the pressure coefficient distributions, but these are really minor changes 
considering the extreme variations in flow topology just described. It should be 
noted, however, that despite being a somewhat unexpected behavior at first, this 
actually mak es some sense if one reasons in terms of boundary layer theory where 
the streamwise pressure variation is determined by the outside flow. It is true that 
this is correct only for unseparated flow, but it can still help the interpretation 
the results obtained. Furthermore, the same basic problems observed in the pres- 
sure coefficient distributions, and discussed earlier in this section, are still present 
regardless of the turbulence model parameters. 

On the subject of obtaining a better pressure coefficient distribution on the body, 
another numerical experiment was performed. Since the data used for comparison 
of the steady state aerodynamic data was obtained for a 7% scale model*- 75 * , there 
was the possibility that the actual flow on the experiment started lamina r over the 
body. The experimental model had carborundum gnt No. 30 applied in a circum- 
ferential band over some portion of the ellipsoidal nose region to ensure transition. 
An attempt was made to reproduce this situation computationally by letting the 
boundary layer start laminar at the nose, and then enforcing the transition at the 
same axial stations where the experiment had the carborundum grit. 

The results obtained indicated that the laminar flow expends faster, and so 
better agreement with the experiment was obtained over the upstream portion of 
the nose region. However, the agreement was already very good over that portion 
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of the body (see, for instance, Figure 6.4). Since the experimental grit was located 
quite far forward on the nose, this beneficial effect of a faster expansion of the 
l»min»r boundary layer was not being felt at the location of the shock, which is 
where the problem was. This, however, suggested the idea of letting the boundary 
layer stay laminar over a larger portion of the nose, perhaps even up to the shock 
location or the point of separation. 

Unfortunately, again, the results were not very encouraging. Although it is 
true that the assumption of laminar flow causes a faster expansion of the flow 
going over the body nose, it is also true that a laminar boundary layer separates 
sooner than a turbulent one. The net result was that the computational pressure 
coefficients would follow the experimental values better over the ellipsoidal nose, 
but the negative peak in Cp was lower because separation would happen sooner. 
In other words, besides being a somewhat arbitrary procedure, this approach was 
causing the shock location, and the separation point, to move too far forward. 

In light of all these observations, the idea was abandoned and the results pre- 
sented in the previous section are for a fully turbulent flow, in the sense that the 
boundary layer is assumed to be turbulent since the very nose of the body. It should 
be mentioned, however, that it is very easy to implement these kinds of tests in the 
current code, since the turbulence modelling routines are probably the only portion 
of the code that is highly modular. Finally, it should be mentioned also that it is 
important for these kinds of calculation to ensure that the flow becomes turbulent 
at some point, at least after separation occurs. Otherwise, the separated flow usu- 
ally becomes aerodynamically unsteady, if we insist in keeping the computational 
solution laminar. It is also unrealistic to expect that separated flow at those flight 
Reynolds number should remain lam i n a r . 
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6.4 Configuration at Angle of Attack: Aeroelas- 
tic Results 

Using the steady state solution described in Section 6.2 as the initial aerodynamic 
solution, aeroelastic analyses were performed for the Atlas-Able IV configuration 
by varying the dynamic pressure (go) while keeping other parameters constant. 
Pure aerodynamic solutions can be performed only in terms of nondimensional 
parameters, if one is interested only in nondimensional coefficients or quantities as 
result. However, when performing aeroelastic analyses, information is necessary 
regarding the (dimensional) freestream values and a reference length in order to 
consistently nondimensionalize the structural dynamic equations. For the same 
reason, it is important to decide whether one is considering a model or the full scale 
vehicle. In our case, we will consider a full scale vehicle since the natural frequency 
information we have in this case is for a full scale model. 

Although the geometrical information about this configuration was obtained 
only in nondimensional form from References [4] and * , the data of Reference ^ ^ ^ 

was dimensional. Considering that the latter presented data for a 7% scale model 
and we wanted dimensions for a full scale vehicle, we chose the reference length 
£o = 0.81280 m , which in this case is the diameter of the afterbody cylinder. 
The freestream speed of sound was chosen to match the experimental data of Ref- 
erence^ 8 ^ and it was set at a n = 336.75 m/s . The freestream density was 
calculated at every run in order to produce the desired freestream dynamic pres- 
sure. Note that by keeping the freestream Mach number and a,, constants, the only 
parameter left to vary the dynamic pressure was the freestream density. 

Also note that, although the speed of sound value being used would corre- 
spond to flight at approximately 1000 m altitude if we consider a standard atmo- 
sphere^^ , the idea cf having a constant speed of sound with a varying freestream 
density is realistic if we think in terms of flight in the stratosphere. In this portion 
of the atmosphere, the temperature is constant and therefore the speed of sound 
is constant, while the density varies with the altitude due to the changes in the 
pressure. It should be also pointed out that we do not have to keep the speed of 

* Ericsson, L.E., personal communication, April 1987. 
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sound constant in the analysis, although we could not devise a more realistic way of 
varying the dynamic pressure while keeping the Mach number constant than the one 
described above. The reason for keeping the Mach number constant is a practical 
rather than a theoretical one. Essentially we want to avoid having to regenerate the 
initial aerodynamic solution for every aeroelastic run, because of the computational 
costs of so doing. The solution is to perform aeroelastic runs varying some other 
parameter, for instance the dynamic pressure, while keeping Mach number, angle 
of attack, and all the other parameters used to generate the initial aerodynamic 
solution constants. 

Three elastic structural mode shapes were used in the present analysis, and the 
rigid body degrees of freedom were assumed to be constrained. The first mode nat- 
ural frequency was estimated as 5 Hz , according to discussions with Ericsson 
A comparison of the first and second mode responses presented by Woods and 
Ericsson 141 , together with some typical launch vehicle responses available from 
Bisplinghoff and Ashley l 60 * , permit the estimation of a value of 17 Hz for the 
second mode frequency. Finally, the third mode frequency was picked at 29 Hz . 

The mode shapes were obtained by an extrapolation of the information avail- 
able in Woods and Ericsson t4] , together with additional information regarding 
the location of the first nodal point for the first mode . The mode shapes 
used are shown in Figure 6.16 . It is important to point out that all modes go to 
zero amplitude at the downstream boundary because this was the way chosen here 
to constr ain the rigid body degrees of freedom. In other words, the downstream 
computational boundary is held fixed throughout the calculation. It must also be 
pointed out that, for xfd greater than approximately 12.5, the mode shapes were 
obtained by extrapolation of the existing data using low order polynomials. More 
specifically, cubic polynomials were used for all three modes, and the boundary con- 
ditions enforced were that they should give the same values of deflection amplitude 
at the last two available experimental points, as well as giving zero deflection and 
slope at the downstream boundary. 

A co nstan t value of rp**» per unit of length was used in these calculations, which 
may be sort of an oversimplification. We do recommend a more careful description 

* Ericsson, L.E., personal communication, April 1987. 
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x/d 

Figure 6 . 16 : Structural mode shapes used for aeroelastic analysis of the Atlas-Able 
IV configuration. 
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of the properties of the vehicle if more accurate results are sought. This 

parameter was obtained from data available on some Atlas boosters , and a value 
of 2555 kg/m was used. The structural damping coefficients were estimated from 
data presented by Woods and Ericsson , and the values used were Ci = 0.0010 , 

C 3 = 0.0018 and < 3 = 0.0036 , for the first, second and third modes, respectively. 

The initial perturbation necessary to start the oscillation is provided, as in the 
cases studied in the previous chapter, by the fact that the initial aerodynamic 
solution is calculated for a rigid vehicle. Since the angle of attack is different from 
zero, there is a net aerodynamic force distribution that is not balanced yet by elastic 
forces due to vehicle deformation. In other words, the body deflected position is 
not in equilibrium with the applied force distribution. 

A typical vehicle response in this case, for an intermediate value of dynamic 
pressure, q D = 400 psf , can be seen in Figure 6.17 , where the time history 
of the generalized modal deflections is shown for all three modes. As it could be 
expected, the first mode response dominates the other modes, and so details of the 
response on the second and third modes are shown in Figure 6.18 . All three modes 
are damped in this case, which is not surprising since we did assume some small 
amount of structural damping and the dynamic pressure is still in the intermediate 
range. What is somewhat unexpected is the fact that the damping coefficient is 
higher t han the plain structural damping coefficient for all the modes, indicating 
that the airstream is contributing positive damping for all the modes in this case. 

This increase in damping is really not very significant for the second and third 
modes, being only of the order of 6% and 3% , respectively. However, the damping 
coefficient for the first mode is 0.0016 , which implies a 60% increase over the 
corresponding structural damping coefficient. It should be noted, however, that we 
are presenting some very small numbers, ic may be questionable how far one can 
trust the precision of the numerical results, considering the amount of computation 
and possible truncation errors involved in these calculations. Nevertheless, the 
qualitative result leaves no margin to question. It indicates that the first mode is 
damped much faster than it would be the case under the structural damping alone. 

It is instructive to mention that the damping coefficient is calculated from the 
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time histories of the response in each mode by assuming an exponential decay. The 
code generated for this purpose initially estimates the mean value of the oscillation. 
Then, using a least squares curve fitting technique, it finds the parameters for the 
exponential envelope of the response. The slope of the natural logarithm of this 
envelope curve is the product of the damping coefficient by the natural frequency 
of that particular mode. From the parameters of the exponetial envelope, the 
structural damping coefficient can be obtained. 

It can also be observed from the results that the frequency of the response in 
each mode is very close to its natural frequency. This is in agreement with the re- 
sults obtained in the previous chapter, where it was observed that the frequency of 
the response remains approximately unchanged regardless of the dynamic pressure. 
This may not seem typical for those used to aeroelastic analysis of lifting surfaces, 
where the mechanism of flutter is usually associated with the merging of the fre- 
quencies for two different modes. However, if an instability exists in the present 
case, it will probably be an one degree of freedom phenomenon, and the fact that 
the frequencies remain unchanged is indeed the expected result. 

The response in terms of the generalized deflections for a much higher value 
of dynamic pressure, qp — 1174 psf , is shown in Figure 6.19 , where all three 
modal amplitudes are presented. As expected, the magnitude of the response in 
the first mode is again much larger than for the other two modes. Details of the 
response in the second and third modes are shown in Figure 6.20 . Except for the 
fact that the amplitudes are about one order of magnitude higher, these results 
are strikingly similar to those obtained for 400 psf dynamic pressure. This is 
not the anticipated behavior. For such high value of dynamic pressure, one would 
expect that the aeroelastic instability condition would be evident, assuming that 
the configuration would develop any aeroelastic problem at this flight condition. 
Even more interesting, however, is the fact that the damping coefficients for all 
three modes in this case are higher than their corresponding values at 400 psf . 

This last observation raised questions about what kind of generalized aerody- 
namic forces were being computed by the method, since a dynamic pressure of 
1174 psf is already somewhat higher than what is actually encountered at this 
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Figure 6.19: Modal response for Atlas-Able IV configuration at 1174 psf dynamic 
pressure (Moo * 0.85 , a * 6*). 
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( b ) 3rd mode. 


Figure 6.20: Detail of the response on the second and third modes for 1174 psf 
freest ream dynamic pressure. 
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Figure 6.21: Generalized aerodynamic forces on Atlas- Able IV configuration at 
Moo = 0.85 , a = 6° and qo = 1174 psf . 

flight Mach number in practice. A plot of the generalized aerodynamic forces for 
the 1174 psf dynamic pressure case can be seen in Figure 6.21 . If we remember 
a previous discussion concerning how the initial perturbation is introduced to start 
the oscillation, the response observed in Figure 6.21 is indeed the expected result. 
In other words, there is a initial transient that takes into account the fact that the 
deflected position is not in equilibrium with the applied aerodynamic load. This is 
shown in the figure by the ini tial drop on the value of the modal generalized forces. 
After approximately 0.05 sec the forces settle down to some nonzero mean value 
together with some small oscillations around this mean. Although not completely 
apparent from Figure 6.21 , this behavior can easily be seen if we look at expanded 
views of the response in each mode after the initial transient has passed. This is 
done in Figure 6.22 . 
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(a) 1st mode. 
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Since the magnitude of the response on the fiat mode is much larger than on 
the other modes, the overall body deflection is dominated by the first mode. There 
is no surprise that the generalised aerodynamic forces follow the first mode oscil- 
lation, with just some minor higher frequency influence. However, the important 
observation is that the force, show very little lag with respect to the deflection. Tins 
can explain why a damped behavior is being obtained even at such high dynamic 
pressures. A study of the modal generalised forces for the 400 psf yielded 

s imilar conclusions. 

A few other case, for different values of dynamic pressure were run, and all the 
results are summarised on the root locus plot shown in Figure 8.23 . As before, 
the free parameter of the locus is the dynamic prereure, and the arrow, indicate 
the direction of increasing dynamic pressure. It can be noticed that the frequency 
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remains approximately unchanged regardless of the value of dynamic pressure, as 

we have previously mentioned. It is evident that no instability is present, and higher 

values of dynamic pressure just seem to move the roots deeper into the left-hand 
plane. 

It is difficult to completely see the behavior of the solution for smaller values 
of dynamic pressure from the root locus plot because of the small variation in the 
frequency. A better understanding can be obtained from a plot of the damping 
coefficient in each mode as a function of the dynamic pressure, which can be seen in 
Figure 0.24 . The first and third modes show a continuous increase in the damping 
as the dynamic pressure is increased. The behavior of the second mode is more 
interesting, showing an increase in damping in the lower range of dynamic pressures, 
then the damping starts to decrease as this parameter is further increased. However, 


115 



CHAPTER «. ANALYSIS OF AN ATLAS-ABLE IV CONFIGURATION 


at the very high value, of dynamic pressure, the damping on the second mode tarts 

to increase again. , , 

From a physical standpoint, the result, indicate that the first and third modes 

show no tendency of presenting aeroelastic problems. The second mode, alt oug 

it never goes unstable, indicates a tendency toward, instability m the dynamic 

pressure mnge of approximately 500 to 800 prf . ft » *•— «»« «“• 

to the result, presented by Woods mid Ericsson 1*1 , the second mode , . the on 

that should present flutter problems. Despite the fact that we are no 

the expected instability, our computational result, are showing some of the correct 

trCI part of the foregoing dbcmpancy between prediction and expectation can be 
explained by the structural data used, which involved some approximations due 
to the lack of the actual vehicle data. It is well known that m*. distribute, 
cm, have significant effect, on flutter result, and, as we have previously dacussed, 
some ratheHever. -options were introduce here. The difficulties mgettmg 
the correct shock strength with the current aerodynamic flow solver may also ha 
some role to play in enplaning part of the difference in the »ults A. « > 
mentioned before, the shock motion mid the shock-boundmy layer intact « 
expected to be an importmit feature of the mechmnmn that causes <■»-**■* 
The capability of numerically reproducing the instability can be impaired by 

lack of a better resolution of the shock region. 

The computational cost, of the* aeroebstic solution, is not trivml, and a few 

comments about this issue would be importmit. For the gnd si* being us m 
the present inveetigation, each aeroelastic iteration take, approximately 12 CPU 
seconds in a CDC Cyber 205 computer. The calculation of the «hicb -P-* 
for each flight condition considered, involves typically somewhere mound 3600 time 
step, due to the limitation, on time step si* in order to keep the tune-accuracy of 
the solution. Of coune, thisisabo dependent upon the range of natural fmquencim 
considered, and the number of iteration, above is correct for the fluencies u«d 
in this work. Actually, it is determined by the lowest frequency, since w, mus get 
at least a few cycle, of the lowest mode in order to fully ascertain the stability of 
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the configuration. 

The significant outcome is that somethin* of the order of 12 CPU hours of a 
Cyber 205 are necessary in order to produce a vehicle response such as ,h, one 
shown in Figure 6.17 . No,, that this doe. no, count the time required to generate 
the initial aerodynamic solution. Furthennore, this is not the overall computer time 
for one aeroeUstic run either. A. we mention* before, due to the heavy amount 
of I/O required on these aeroelastic solutions (because the datable is not core 
contained) , the CPU time represents only 30% of the overall system time. 


6.5 Study of a Zero Angle of Attack Case 

The Atlas-Able IV configuration was also studied at a zero angle of attach flight 
condition. Since the geometry was the same, we used the srnn. grid system previ- 
ous‘y generated for this analysis. The structural and aerodynamic data were also 
the same, except tha, in this case the angle of attack was set to zero, and so were 
all the structural damping coefficients. The initial aerodynamic solution was, then, 

* ‘ fW * he ~ 0 85 “ d ° = °* ftgkt condition. The Reynolds number 

was the same as before, and the boundary layer was considered turbulent 

A flow visualization of ,h. steady state solution obtained, which will be the 
starting solution for aeroelastic analysis, can be .sen in Figure 6.25 . The figure 

IT M “ h nU “ Ur “ d Plot, for a side view 

of the body. The solution is asymmetric, as is clear from the figure since th, 

top snd bottom contoun are a mirror image of each other. Even a, zero angle of 

“ " fl ° W “ Pini “ on ow «“ This is evidenced by th, reversed 

e oci y pro es shown in Figure 6.26 for a portion of a longitudinal plane. The flow 
separate still over the forebody cylindrical section, and the boattail is completely 

unmer«d m the region of reversed flow. Reattachment will occur only downstream 
of the boattail over the afterbody cylinder region. 

The computed pressure coefficient distribution over the vehicle is compared to 
expennumtal result, 51 as a way of validating th, ,te*y .tote nondynamic solu- 
■on. Tins companson can be seen in Figure 6.27 . I, is no, particular important 
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Figure 6.25: Flow solution about an At 
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Figure 6.27: Pressure coefficient distribution over an Atlas-Able IV payload at 
Moo = 0.85 and a = 0* . 
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which longitudinal plane of results was used for this comparison, because the solu- 
tion is axisymmetric. However, for the sake of completeness, the plane correspond- 
ing to the leeside in case of a positive angle of attack was the one used. The results 
show the same features observed in the solution at angle of attack. The compu- 
tational pressure coefficients foUow the correct trend of the experimental results, 
but the shock resolution is less than perfect. As was previously discussed, problems 

with the numerical dissipation algorithm may be causing the discrepancies in these 
results. 


Only one aeroelastic case was run at this flight condition because of time limi- 
tations. The dynamic pressure was set at 400 psf . The structural data used were 
the same considered in the cases studied in the previous section, except that the 
structural damping coefficients were set to zero in the present simulation. Since the 
angle of attack is zero, we need a different way of introducing the initial perturbation 
that starts the oscillation. An initial displacement of the first mode generalized co- 
ordinate, while keeping the other two modes undisturbed, was the method selected 
here. The magnitude of the initial deflection was selected to be of the same order 
of the maximum elastic displacement observed in that mode at 6° angle of attack 
when the same dynamic pressure was being used. 


Tune history of the response in terms of the generalized modal deflections is 
shown in Figure 6.28 . It must be emphasized that all three modes are shown on 
this figure. However, due to the large difference in the magnitude of the response, 
the second and third modes are only seen as a zero line. An expanded view of these 
latter modes response is presented in Figure 6.29 . Just by looking at the results, 
it is clear that the oscillation amplitude of the second and third modes seems to be 


growing. The first mode appears to be neutrally stable. Furthermore, the higher 
modes seem to be modulated by by the first mode response. 


A more complete analysis of the first mode response can be accomplished using 
the same software developed for the study of the previous aeroelastic runs. The 
result of the analysis is that the response is at a frequency of 5 Hz , which is 
exactly the natural frequency for this mode, and that the damping coefficient is 
^ 9 x 10 6 . In a strict sense, this mode would be called unstable since the 
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(6) 3rd mode. 


Figure 6.29: Expanded view of response on the second and third modes at 
Moo = 0.85 , a = 0° and qo = 400 psf . 
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damping coefficient is negative. However, considering the smallness of the value 
of Ci , it is probably more appropriate to accept the mode as neutrally stable. It 
is important to look at the problem from a practical point of view too. Since the 
vehicle will be accelerating through the transonic regime, and not flying continuously 
in it, such an extremely slow-growing instability would probably be of no concern. 

Due to the apparent modulation of the response on the second and third modes, 
their study is a little more complex. We decided to study first the frequency content 
of the responses shown in Figure 6.29 . This can be done by taking the Founer 
transform of those time traces. The result of such operation is shown in Figure 6.30 , 
where we are plotting the magnitude of the FFT of the response on each of those 
modes. The reader should note that, rigorously, the Fourier transform of a diverging 
signal does not exist. Here, however, we are simply taking a discrete FFT over a 
finite time interval. For the second mode response, the peaks at 5 and 17 Hz are 
evident from Figure 6.30(a) . Similarly, for the third mode, peaks at 5 and 29 Hz 
can be seen in Figure 6.30(6) . 

It is interesting to note the relative magnitude of the response in each frequency. 
The response on the third mode seems closer to the expected result, where the 
spike at 29 Hz is much larger than the one at 5 Hz . The trend is reversed in the 
second mode response, where the 5 Hz peak has a larger magnitude than the one at 
17 Hz . This explains why the third mode response is much “cleaner”, in the sense 
that it is much easier to distinguish between the higher frequency response and 
the modulation itself. There seems to be more interaction of the two frequencies 
in the second mode response, which accounts for the more complex pattern in 
Figure 6.29(a) . Moreover, with the two frequencies being closer to each other, it is 
understandable that the time trace should get more complex. 

Although it clear from Figure 6.29 that the amplitude of response in 

those modes is increasing, we attempted to filter out the 5 Hz portion of it. This 
be accomplished by eliminating the 5 Hz component of the Fourier transform 
of the response, and then taking the inverse transform. We were not completely 
successful in eliminating the modulation, and the filtered third mode response does 
not add much information to what can already be seen from Figure 6.29 (6) . For the 
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Figure 6.30: Frequency content of the modal response of an Atlas-Able IV vehicle 
at dynamic pressure qo = 400 psf . 
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Figure 6.31: Filtered second mode response for Atlas- Able IV vehicle. 

second mode, however, the filtered response clearly indicates the growing amplitude 
of oscillation. The latter is shown in Figure 6.31 . 

As we have seen from the results of the previous section, it is also interesting 
to look at the modal generalized forces. These are shown in Figure 6.32 . There is 
no surprise that the forces follow the first mode deflection, since the overall vehicle 
deformation is almost completely deter mined by the first mode. It can be observed 
that there is very little phase lag between the motion and the aerodynamic forces. 
Figure 6.33 shows the aerodynamic load distribution along the body for several 
instants during the oscillation. Since, again, the first mode amplitude dominates 
the other ones, we used the first mode generalized coordinate as the total deflection 
amplitude when selecting at which instant of time we would sample the running 
normal force. Essentially, Figure 6.33 is showing that most of the load is located 
at the vehicle elliptical nose for the zero angle of attack case. This explains why 
the magnitudes of the generalized forces are so similar, since the mode shapes are 
similar to one another for that body region (see, for instance, Figure 6.16). 
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Figure 6.32: Modal generalized forces on Atlas-Able IV at M* = 0.85 , a = 0° and 
q D = 400 pef . 
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Figure 6.33: Unsteady aerodynamic load distribution on Atlas-Able IV at various 
points along the first mode oscillation. 
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Finally, as we mentioned before, the second mode was expected to become aero- 
elastically unstable, according to results presented on Reference (4] . For the zero 
angle of attack case, our computational results are capturing the instability appro- 
priately, as evidenced by the foregoing analysis. It is also interesting to note that 
our predictions regarding the proneness of the smaller angle of attack situations 
to aeroelastic instability were confirmed by the results. Despite the fact that our 
calculations at the 6° angle of attack condition assumed nonzero structural damp- 
ing coefficients, we saw that at 400 psf dynamic pressure the airstream was still 
contributing positive damping. On the other hand, an unstable solution is obtained 
at zero angle of attack for the same dynamic pressure. So, at least qualitatively, 
our computations are predicting the correct behavior in the cases analyzed. 
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Conclusions 


7.1 Summary 

The idea of using computational fluid dynamics techniques for aeroelastic analysis 
in the transonic regime of flight is not new. Several of the references cited through- 
out this text substantiate this statement. In this context, the contnbutions of the 

present work are associated with the use of the three-dimensional unsteady 

pressible Navier-Stokes equations for the flow simulation and aeroelastic analysis of 
realistic 3-D body configurations. It is a substantial departure from past investiga- 
tions, where a potential formulation was commonly used to solve the aerodynamic 
problem or where only two-dimensional flows were considered. 

The present research initially studied a hemisphere-cylinder configuration. Steady 
state, as well as unsteady aerodynamic, results for this geometry were presented 
and di iscuased. The steady state results showed good agreement either with exist- 
ing experimental information or with data obtained from other comparable CFD 
codes. They provided good insight into the code’s behavior for both subsonic and 
supersonic freestream conditions. No aeroelastic cases were investigated for the 
hemisphere-cylinder. On the other hand calculations for forced rigid body oscilla- 
tions were performed, wi u ihe primary objective of acquiring a better understand- 
ing of the parameters involved in these unsteady computations. 

The unsteady airloads due to this forced, sinusoidal rigid body oscillation could 
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not be compared to any results in the literature because we were unable to locate 
similar experimental investigations. To the author’s knowledge, this was the first 
attempt to capture the unsteady behavior of such a configuration undergoing sinu- 
soidal pitching oscillations using a Navier-Stokes formulation. Happily, the unsteady 
calculations reproduced the expected flow features. Due to time constraints, how- 
ever, we refrained from studying cases at lower reduced frequencies which probably 
would be more representative of actual flight conditions. 

The study of a general hamm erhead payload provided further confidence in the 
code developed by producing physically reasonable results for a transonic flow con- 
dition. Hammerhead vehicles are by no means simple geometries for computational 
simulation, because some form of flow separation is usually associated with the 
presence of a boattail in the forebody. This is in addition to the shock-induced sep- 
aration due to the transonic shock impingement on the body surface. The steady 
state aerodynamic results reproduced the expected flow features. The aeroelastic 
calculations showed that the first configuration studied was free from instability at 
the flight conditions considered. 

The stable aeroelastic behavior of this vehicle can be explained by the exis- 
tence of a longer hammerhead section, i. e. , a longer cylindrical forebody region. 
It prevents the merging of the two separation regions mentioned above, therefore 
resulting in a milder unsteady aerodynamic environment This observation agrees 
with existing empirical rules for the design of hammerhead payloads. Besides the 
physical considerations of aeroelastic stability of this configuration, its study was 
instrumental in developing the details of the structural-dynamic portion of the code. 
A more appropriate assessment of the computational requirements for aeroelastic 
analysis was possible, and acceptable time step levels could be selected. 

For the Atlas-Able IV vehicle, steady state wind tunnel pressure coefficient dis- 
tributions were available in the literature for various t unn el conditions in the tran- 
sonic regime. Actual in flight transonic aeroelastic stability information was also 
available, as well as the results of semi -empirical analyses. The details of the exact 
flight conditions encountered were, unfortunately, rather sketchy. The configura- 
tion was studied at two different angles of attack, and for each flight condition the 
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dynamic pressure was the parameter varied in the search for the stability boundary . 

The steady state results compared well with experiment, except that a better 
computational resolution of the shock region would be desirable. In the process 
of perfor ming these calculations, some numerical difficulties with the present flow 
solver were uncovered. The artificial dissipation algorithm used showed some lim- 
itations in the sense that, in order to maintain numerical stability, the required 
magnitude of this parameter was already interfering with the physics of the flow. 
The turbulence model employed proved very sensitive to some of its parameters, 
producing a wide variety of flow topologies for the same flight condition when the 
only thing varied was the distance of the search for the maxim um velocity in the 
boundary layer profile. 

This last observation can probably be explained by the fact that the kind of 
massively separated flows which occur in this case may exceed what the turbulence 
model was originally derived for. A discussion of the possibilities for further upgrad- 
ing the current flow solver algorithm will be presented in the next section. While on 
the subject of flow solutions, it must be said that the availability of a good graphics 
post-processing package is absolutely essential for interpretation and understanding 
of these three dimensional CFD calculations. The study of body pressure distribu- 
tions only, even though important, is not enough to ensure that the correct solution 
is being obtained. A very good example was seen in the numerical experiments 
performed with the turbulence model parameters. 

The aeroelastic calculations reproduced the known Atlas-Able IV results qualita- 
tively very well. The second mode instability ^ was clearly observed at zero angle of 
attack and M* = 0.85 . At an angle of attack of 6 degrees and the same freest ream 
Mach number', no unstable aeroelastic behavior was predicted. However, the second 
mode again showed some tendency of becoming unstable for some portion of the 
dynamic pressure range investigated. The instability, when it occurs, is certainly a 
one degree of freedom phenomenon. The frequency of the response in each mode 
remains approximately unchang ed, regardless of the variations in dynamic pressure. 

The Atlas-Able IV configuration has a much shorter cylindrical forebody sec- 
tion than the first hammerhead configuration analyzed, which explains the merged 
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separation. As previously discussed, in these cases the shock-induced separated re- 
gion merges with the flow separation due to the adverse pressure gradients at the 
boat tail. The result, in terms of vehicle stability, is that the possibility is created 
of larger lags in the airloads with respect to the motion. Therefore, the risk of 
aeroelastic instability is increased. As extensively discussed in the previous chap- 
ter, the computational results also confirmed the prediction that the smaller angle 
of attack cases would be more susceptible to instability. However, despite the good 
qualitative agreement, improvements in the aerodynamic solver and more accurate 
structural data are still necessary for a quantitative prediction of flutter speeds by 
means of the present method. 

The CPU time per iteration for aeroelastic analysis is only 20% higher than 
for steady state aerodynamic calculations. The time spent on purely structural- 
dynamic calculations in the present implementation is, however, truly negligible. All 
the computational time increase is associated with the need to regenerate the grid 
at every time step when performing aeroelastic calculations. Typically, approaches 
based on modal analysis will spend almost negligible time on the solution of the 
structural-dynamic equations. On the other hand, if the vehicle is represented by 
some finite element model, this time is certainly bound to increase. 

It must be emphasized that the increase in overall system time per iteration 
is much higher for aeroelastic than for purely aerodynamic examples. Since the 
database is not core contained, some rather extensive I/O is required every iteration 
even m steady state calculations. Use of asynchronous I/O allows an efficiency of 
60% for pure aerodynamic computations. Efficiency here is being measured as the 
ratio of CPU time to overall system time. For aeroelastic analyses, however, the 
code efficiency is reduced to only 30% , because of the considerable increase in 
the input/output required per iteration. These efficiency figures are, of course, 
very machine dependent. The discussion here applies only for the Cyber 205 [ 

although other machines may experience similar trends for databases that are not 
core-contained. 

It is true that the computational requirements for using of the present method 
are still quite high. However, considering recent improvements in hardware as well 
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as numerical methods, it is fair to say that we are getting very close to acceptable 
levels. On the the storage side, i. e. , in terms of the internal memory required, 
environments such as the Cray 2 supercomputer completely eliminate the need for 
out-of-core temporary storage. Our requirements, which are for about 7 million 
words in the database plus a comparable amount for temporary arrays, can be 

considered small for a 256 mega-words core machine. 

Since in most systems a problem of this size would require reads/writes to tem- 
porary files, the subject of computational time must be looked at on the two levels 
of CPU time and overall system time. There is no question that, as far as the latter 
is concerned, the optimum condition is obtained when everything is kept in core . 
The problem with I/O sometimes can be minimized by intelligent structuring of the 
database and by streaming the flow of information in the code to take advantage 
of that. As the statistics given in the previous chapters show, the present code was 
successful in doing this for steady state problems. For aeroelastic analyses, however, 
the performance dropped considerably. 

A reduction in actual CPU time can be accompliahed either by faater flow solver 
algorithms or by faster machines. The search for faster and more reliable computa- 
tional schemes is a continuing objective of CFD research. A breakthrough in terms 
of raw machine speed was realised with vectorisation. More recently, some attention 
has been given to the approach of having a multi-processor* computer, where each 
processor could work on a portion of the program at the same time. An example 
of this technology is the ETA-10 system. It is worth pointing out that the pencil 
data structure used in the present work is very suitable for such multi-processor 

environments. 

A discussion of possible extensions for the method is now presented. The main 
purpose is to relate the work done here with similar applications, which may in- 
volve even more complex configurations or flight conditions than the ones described 
in the preceding chapters. The aerodynamic formulation is very general. Actu- 
ally, it is as general as one would require for the flight regime and vehicle types 

1 Note that the nomenclature may be misleading, since there are many machines which have m 
than one processor but do not have the capability we are describing here. 
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under consideration, except with respect to the representation of turbulence. The 
structural-dynamic formulation, however, has several assumptions which may seem 
quite restrictive at first. These assumptions are believed adequate for the physical 
situation being treated by the present work. Nevertheless, a user may be interested 
in generalizing this formulation. 

If we assume that a modal approach would still be used, the major general- 
izations that could be incorporated would be to introduce rigid body degrees of 
freedom, to consider that the modes might be elastically or inertially coupled, and 
to allow for a full three dimensional beam type of behavior. The present theoretical 
formulation already allows for rigid body degrees of freedom. However, they were 
never implemented in the code, and a few programming changes would be required 
for their introduction. 

If normal vibration modes are not known, i. e. , the modes to be used are 
coupled either elastically or inertially, the method can be easily extended. The 
difference would be that, instead of solving for each of the generalized coordinates 
separately, the solution would have to be done considering the (now) coupled system 
of equations. In terms of computational effort, this means that matrices would have 
to be inverted, as opposed to the solution of a number of scalar equations. One 
should note, however, that the number of structural modes considered in a typical 
aeroelastic analysis is extremely small compared to the size of the matrices involved 
in the solution of the aerodynamic problem. As a result, the overall computational 
cost of the solution should not be affected very m uch . 

The aerodynamic formulation allows for arbitrarily large angles of attack, which 
includes cases where asymmetric separation may occur. It would be interesting to 
give to the structural-dynamic formulation the capability of handling the same flight 
condition. This would involve considering bending in two planes and possibly tor- 
sion of the body about its longitudinal axis. This extension is very straightforward 
in concept, although its programming may require some very careful thought. 

On the other hand, if the current low angle of attack restrictions are kept for the 
structural-dynamic formulation, some simplifications can be incorporated in the flow 
solver. The computational grid used in all the calculations presented here is a 360° 
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mesh in the circumferential direction. The results presented in the earlier chapters 
show that the flow solutions obtained are symmetric with respect to the pitching 
plane. In other words, the pitching plane can be used as a plane of symmetry for the 
flow solution, and almost half of the grid points could be saved. This simplification 
would cut the complete database in half, and the savings in computational time 
would be much larger since this time increases nonlinearly with the size of matrices. 

The use of an algebraic grid generation scheme proved to be appropriate for the 
applications in the present work. However, this may not be true if --ore complex 
geometries are considered. There are some hammerhead configurations with very 
steep boat tail angles. In such cases, the approximation used here of letting the 
normal grid lines come into the body surface at angles slightly different from 90° 
may be inadequate. Although algebraic schemes could still be used for these cases, 
more care would have to be exercised to ensure that normal grid lines do come 
in normal to the body. Some launch vehicle configurations have fins, or strap-on 
boosters, or even cluster designs. For these cases, certainly more elaborate grid 
generation schemes would be required. 

The procedure for locating the downstream boundary followed in the present 
approach was possible because all the physical phenomena of interest were hap- 
pening at the forebody. If configurations like the ones mentioned above are under 
consideration, the loads on the afterbody may significantly influence the overall ve- 
hicle’s aeroelastic behavior. The more appropriate strategy for those cases would 
be to introduce the body base into the computational domain. The downstream 
computational boundary would have to be moved a least one body length down- 
stream. This would make the flow solution more complex, in the sense that mixing 
layers and base flows would appear in this downstream portion of the computa- 
tional However, since we are not interested in these phenomena as far as 

our aeroelastic analysis is concerned, they would have to be resolved only to the 
extent of avoiding significant inaccuracies in the aerodynamic loads on the body. 

In conclusion, the method developed in the present work seems to be an at- 
tractive option for aeroelastic analyses of launch vehicles in the transonic phase of 
flight. In this regime, the usual linearized methods break down, and an approach 
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that would take into consideration the aerodynamic nonlinearitiee is required. The 
applications presented provide confidence in the method, but it is clear that some 
details still need improvement before the code could be used in a production envi- 
ronment. There are launch vehicle, being currently demgned (e. g. , a commercial 
tlaa/Centaur and the Titan IV series) which probably could benefit from the tech- 
mques presented here in terms of reducing design coats. 

It should be kept in mind that in the recent yean there ha, been a subetantial 
reduction in the cost of computation, whereas experiments are still very expensive, 
t is the author’s belief that, eventually, computational procedures like the one 
proposed here will be capable of full integration into the design process. Thus, 
t ey may permit a more rational aeroelastic design instead of the usual fixes that 
characterized the aeroelastic clearance of flight vehicles over the years. 


7.2 Recommendations for Future Work 

The applications presented throughout the work exemplify the capabilities of the 
method. Despite the good qimlitative agreement of the prment results with existing 
data. It was clear that further consideration of. few computational problems was 
s .11 necessary. The difficulties encountered with the artificud dissipation model 

d ,T™ OTmi “ liOT - F ” » comparison of the magnitudes 

° J V “ C08lty COefficiCnt of the amount of artificial dissipation being 

mtrod W would be instructive. Moreover, the question left answered in the present 

wor o w y we are being forced to use such large amounts of artificial dissipation 
m order to keep the scheme stable should be resolved. 

* “ IUtion 10 ,his P roblem “Sfct "imply be that a more elaborate 

artificial dissipation scheme is necessary. Recent versions of the ARC3D code! 52 ! 

have a nonlinear numerical dimipation model The amount of artificial dissipation 
introduced at each computational point is weighted by a spectral mdiu. soding 
whKh is a sum of the spectral mdii of the Jacobian matrices. This model produced 
better results, and it may be helpful for solving the difficulties encountered by our 
present code with respect to artificial dissipation. 
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The foregoing discussion assumes, however, that the Beam and Warming ap- 
proximate factorization scheme! 57 ! will continue to be used. This does not have to 
be true. Recent so-called upwind scheme, eliminate the need to explicitly introduce 
artificial dissipation into the scheme. The numerical dissipation is already built 
into the scheme by the way the governing equations are approximated by means of 
one-sided difference formula The replacement of the current solver by these up- 
wind schemes may be a possible form of improtring the present code. Furthermore, 
so-called TVD schemes 179 ' are well known for having very good shock capturing 
features. A possible drawback is that these schemes may be slower than the Beam 


and Warming algorithm. 

As we have seen, another point of serious concern was the turbulence mode 
Some work could be done in this area in order to enhance the present flow solver 
algorithm. The tradeoff used to be that more elaborate models, despite produc- 
ing better results for some flow cases, still would not yield significant enough im- 
provements in most awes to justify the mcre.se in complexity mid computational 
requirements. Recently, however, some turbulence models have been created (e. g. , 
the Johnson and King model' 80 ') which advance the numeriad representat.on of 
turbulence with very little increase in computational coeta. 

Another enhancement to the present code could be accomplished by tncorpo- 
rating some self-adaptive grid techniques' 81 ’ 821 into the current algorithm. These 
techniques let the position of the grid points be driven by the solution, concentrat- 
ing them on region, of higher flow gradients. The tradeoff usiwlly involved is that, 
by portioning the grid point, in some optimum way, fewer point, are oecewwry for 
a given accuracy. Thus, a faster run time is achieved. On the other hand, tune is 
consumed to si.pt the grid, which increwes the total run time. Since in the present 
approach the grid is being regenerated at every time step for aeroelas.ic anTysw, 
we would only get the advantages of self-adaptive grid method, without paying the 

additional costs. 

The subject of improved graphic capabilities for vumalization of flow results has 
been extensively discussed in the literature. The author believes that the current y 
available packages («. g. , PL0T3D, RIP, GAS, etc. ) are very useful for the 
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analysis of steady state aerodynamic results. When dealing with unsteady problems, 
however, the requirements are more stringent. Static colored figures are not enough, 
and one has to resort to animation techniques in order to comprehend fully all the 
phenomena involved. The amount of data that has to be stored and then shipped 

across the network in order to use the existing graphic and animation packages is 
tremendously large. 

As a result, there is a tendency to settle for an analysis of modal generalized 
deflections, as we did here, or some integrated aerodynamic quantity. Although 
these are important too, there is much more information being generated by the 
computation that is never studied because of the lack of practical ways of displaying 
it. For instance, a visualization of how the flow topology is modified as the body 
oscillates might be very helpful for understanding the mechanism behind a possible 
aeroelastic instability. Although the necessary tools for these visualizations exist 

today, their use is still so time consuming as to discourage a more widespread 
utilization. 

The current trends in launch capabilities have given rise to a renewed interest 
m expendable launch vehicles. As discussed above, hammerhead payload config- 
urations will certainly be under consideration. This fact should provide enough 
motivation for the continuing development of aeroelastic analysis techniques like 
the one presented here. Furthermore, it should also provide plenty of experimental 

data that could be used to validate these techniques and to guide their develop- 
ment. 
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Jacobian Matrices 


The three-dimensional, inviscid flux Jacobian matrices in generalized curvi lin ea r 
coordinates are given by 

K t ** 

Kg<j> 2 — uB Kf + 6 — Kg (“7 — 2) U 

Ky<j> 2 - vB K s v - «„ (7 - 1 ) U 

k, 4 2 - wB KgW - k m (7 - 1) U 

B (24? - yep' 1 ) Kg {yep' 1 -<?)-{ y-l)uO 

Kg 0 

KyU — Kg(y — l)v Kgti — Kg{ 7 — l)tu K x (y — 1 ) 

K t + B - Ky(y - 2)v K t v - Ky{y- l)ti> (A- 1 ) 

KyW — fC, (7 “ 1) v *1 + B — K* (7 “ 2) tl> K x { 7 — 1) 

Ky{yep- X -<j> 2 )-{y-\)vB K.iycp ' 1 - <f> 2 ) - {y - l)wB K t + yB . 

where 

B = «*u + « v t> + *,to 

= |(7 “!)(«* + »* + »*) 

Here * = {, jj or C for A, B or C, respectively. 
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The viscous flux Jacobian matrices are obtained from 

0 0 0 0 o 

<*iMp -1 ) <*2 6 „(p _1 ) a i S H (p- 1 ) 0 

”■*31 q 2 a A 6 K {p-') a h 6 H {p- 1 ) 0 (A. 2) 

m « ai6«(p-') <**6. (p~ l ) 0 

m 5i 771 S2 m 53 rr»54 a 0 S lt (p~ 1 ) 

where 

m 2 i = ~ a i£« ( u /p) — ojS H (v/p) — ar 3 6 „ (w/p) 

m 31 = ~ a 2 ^« (u/p) — Oli^H (u/p) — a$S H (t u/p) 
m <! = _a 3 ^« (w/p) - ar 5 <5„ (v/p) - ar«£„ (w/p) 
m 51 = -ao6 H [(e/p 3 ) - (u 3 + u 2 + u 2 ) /p] 

= (w 2 /p) - a 4 <5„ (u 3 /p) - a«£„ (u? 2 /p) 

= - 2 a 2 £„ (uv/p) - 2 a 3 £„ ( uw/p ) - 2 a 5 £« ( vw/p ) 
m 32 = -Oo£« (u/p) - m 21 
”353 = -Oto6* (u/p) - m 31 

m 54 = -a 0 <$„ (u/p) - m 41 
a 0 = 7 (Pr * 1 + p t Pr~ l ) (* 2 + K * + * 2 ) 

a » = ( 1 +P«)(|^ + «J + ^) 

1 „ 

°2 = g(l+P«)*r*„ 

^3 38 ^(1 fit) 

<** - (i + p t ) 

a * = jt 1 + Pt) K V K * 
a « = (1 + £<) (*2 + *; + |* 3 ) 

As in the inviscid case, « a £, tj or C for A/ << A/„ or Af <t respectively. Furthermore, 
we remind the reader that all variables in these expressions are dimensionless. In 
particular, p< Pt/p oo> where ^i w is the freest ream (laminar) viscosity coefficient. 
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Periodic Tridiagonal Solver 


A periodic, block tridiagonal matrix has the general form 

B\ Cj 0 • • • 0 Ao 

B 7 C 3 ••• 0 0 

0 Aj fl3 '■ 

0 ’ • C»-j 0 

0 : Bn—\ Cn 

C n+ 1 0 0 A n -1 


(B.l) 


For the present appUcation, each element in A is a 5 x 5 matrix. However, there 
is another “dimension” to A, which we are not showing here. Since the domain 
of calculation is a parallelepiped in computational space, we actually operate in 
a whole plane of data when performing the solution in a certain direction. This 
other dimension of A, which one may call its “depth”, is equal to the number 
of grid points in this plane. The reader should also realize that we are using a 
somewhat loose nomenclature by calling these entities “matrices”. They would be 
(more appropriately) referred to as multi-dimensional arrays. 

It is instructive to think in terms of matrices because they are easier to represent. 
Furthermore, for each point in this grid plane (or grid surface, if one thinks in terms 
of the physical space), the finite difference equations are matrix relations. On the 
other hand, it is important to remember that we are solving for a whole data plane 
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at once, since the storage requirements are determined by the size of the resulting 
multi-dimensional arrays. 

The process of inversion involves an L-U decomposition of matrix A, such that 


A = CU 


where 

L\ 0 0 • • • 0 0 

ij 0 • • • 0 0 

0 M % Li •. i - 
: : 00 
0 0 i •. L n . x 0 

t Ni ••• Nn-2 M n -x L n 

I U-i 0 ••• 0 Vi 

0 / U 3 ••• 0 V 7 

0 0 1 ••• 1 

0 0 0 U n -x v n _ 7 

: Un 
00 0 • •• 0 1 

Here I is the identity matrix. 


(B.2) 


(B.3) 


(B.4) 


In order to present the algorithm for solution adopted in the present work, we 
will assume that the solution for a system of the form 


AX = Y 

is being sought. Furthermore, we will call the intermediate result 

X = UX = C~ l Y 


where 



(B.5) 
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The solution algorithm can be written as follows: 


Forward sweet 


• First element (i = 1): 


L x = Bx 
Mx = Ax 
Nx = Cn+l 

lx = L?yx 
Vi = L?Ao 

T. = NxVx 


• For i = 2 until n — 2 : 


Ui = L'lxCi 
L { = Bi-Mi-xUi 
Mi = Ai 
Ni = -Ni-xUi 
7i = I“ l (yi - Mi-ili-i) 
V; = -L-'Mi-xV-x 

t. = t. + aw 


• Next-to-last element (* = n - 1): 

Ui = L~lxCi 
Li = Bi-Mi-xUi 
Mi = Ai-Ni-xUi 

Ti = ir 1 (y< - Afi-iTi-i) 

= Mi-iV;.! 
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• Last element (i = n): 

U n = L-l, (C n - T h ) 

r. = r. + M n . x u n 

Lrx = B n - T a 

~ L n l (y n — Mti-l^n-l) 

Backward sweep 

• Last and next-to-last elements: 

= J»» 

^n—l ~ ^n— 1 U n X n 


• For i = n — 2 until 1 : 


^i+l^i+l V»X n 

As the expressions above show, only the arrays required for the backward sweep have 
to be completely stored. At a given point along the forward sweep, only the arrays 
in the current col umn of C, and i!n the preceding one are needed. All information 
prior to that can be discarded. Furthermore, the arrays in A are generated as they 
are needed during the forward sweep. In summary, care is exercised in order to 
reduce the storage requirements to a minimum. Further reduction in storage space 
could only be achieved at costly penalties to the computational time, since it would 
involve regenerating portions of U during the backward sweep. 
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